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ABSTRACT 

The  problem  studied  in  this  thesis  is  the  guidance  of  interplane- 
tary vehicles  which  are  thrusting  for  a  large  portion  of  the  transfer. 
The  vehicle  is  represented  by  a  seven  component  state  vector  consist- 
ing of  the  position,   velocity  and  mass  of  the  spacecraft.     The  analysis 
is  linearized  by  assuming  that  the  actual  state  of  the  vehicle  differs 
only  a  small  amount  from  a  known  reference  state.     The  reference 
trajectory  is  assumed  to  be  a  propellant- optimal  path  connecting  the 
initial  and  final  points. 

The  goal  of  the  postulated  guidance  system  is  to  satisfy  position 
and  velocity  conditions  at  the  target  with  minimum  propellant  expedi- 
ture.     Both  fixed-time-of-arrival  and  variable-time- of-arrival  guid- 
ance are  discussed.     Specification  of  the  guidance  criterion  in  the 
above  manner  permits  the  techniques  of  optimal  control  theory  to  be 
applied  to  the  problem.     Emphasis  is  placed  on  finding  an  analytic 
solution  of  the  linearized  equations.     The  desired  solution  is  the  con- 
trol program  which  satisfies  boundary  conditions  and  minimizes  pro- 
pellant expenditure. 

The  method  for  solving  the  guidance  problem  is  shown  to  be  suit- 
able as  a  technique  for  computing  optimal  reference  trajectories.     The 
trajectories  are  computed  by  iterative  application  of  the  guidance  sol- 
ution.    Application  of  the  guidance  solution  to  the  trajectory  problem 
is  shown  to  exploit  an  interpretation  of  the  Euler  equations  which  per- 
mits simplification  of  the  computation  technique. 

The  guidance  solution  is  tested  in  a  numerical  example  by  using 
it  to  compute  trajectories  from  Earth's  orbit  to  the  Martian  orbit  for 
different  low-thrust  vehicles. 

The  guidance  solution  is  based  on  the  assumption  that  vehicle 
state  is  known  at  the  time  a  new  control  program  is  to  be  generated. 
Prior  studies  by  several  investigators  detail  methods  of  using  celestial 
measurements  to  estimate  state.     A  portion  of  this  report  is  devoted 
to  extending  the  method  of  celestial  measurements  to  include  measure- 
ment of  engine  performance.     The  additional  measurement  is  shown 
to  improve  the  estimate  of  state. 
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otherwise  noted.  A 
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Example 

2.  Underscored  letters,   both  upper  and  lower 
case,   represent  column  vectors.  P,    r 

3.  Superscript  T  represents  the  transpose  T        T 
of  a  matrix  or  vector                                                            A    ,    P 

4.  Superscript  -1  represents  the  inverse  - 
of  a  matrix                                                                                A 

5.  Unless  underscored,   lower  case  letters 
represent  scalar  quantities  r 

6.  Juxtaposition  of  matrix  and  vector  sym-  A  B  „ 
bols  represents  matrix  multiplication  PQ 

70   The  determinate  of  a  matrix  and  the 
magnitude  of  a  vector  (when  lower  case 
letters  are  ambiguious)  will  be  indicated 
by  vertical  bars.  |A|or|rj 

8.  Vertical  brackets  indicate  a  column 
vector  composed  of  the  enclosed 
quantities. 

9.  Square  brackets  indicate  a  matrix  ,  n 
whose  elements  are  the  quantities  la  b]  orJa  b[ 
enclosed                                                                                    |_c  dj       (c  dj 

10.  Diamond  brackets  indicate  the  time  T 

average  of  the  enclosed  quantities.  <  e    e     > 

The  conventional  dot  notation  is  employed  to  indicate  the  time 
derivative  of  a  quantity  with  respect  to  a  non- rotating  reference  frame. 

Subscripts  are  used  to  supplement  the  fundamental  notation. 
Subscripted  variables  are  defined  as  they  are  introduced. 
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CHAPTER  I 


INTRODUCTION 


Early  in  the  investigation  of  rocket  propulsion  for  extra-terrestrial 

travel,  it  became  apparent  that  conventional  chemical  propellants  were 

12  3  4 
inadequate  for  many  interesting  space  missions  .  This  fact  is 

due  to  the  relatively  low  energy  content  per  unit  mass  (specific  energy) 
of  chemical  fuels.      The  theoretical  mass  ratios  required  for  many  inter- 
esting missions  in  the  solar  system  approach  numbers  of  the  order  of 

3  5 

10     and  higher  when  chemical  fuels  are  used 

The  need  for  more  efficient  means  of  propulsion  has  led  to  the 
investigation  of  energy  sources  other  than  chemical  reactions.     Such 
studies  have  produced  an  entire  spectrum  of  propulsive  techniques    , 
each  having  particular  mission  capabilities  and  each  having  its  own 
theoretical  and  practical  difficulties.      In  general,  the  use  of  higher 
specific  energies  is  accompanied  by  a  decrease  in  propellant  flow  rate 
and  longer  propulsion  time    for  a  given  mass  ratio.     In  limiting  cases 
the  propulsion  time  equals  the  transfer  time.     In  addition,  a  longer  pro- 
pulsion time  permits  reduction  of  thrust  levels  for  a  given  total  impulse. 
The  terms  "continuous  thrust"  and  "low  thrust"  stem  from  these  two 
effects  of  high  specific  energy.     Although  the  terms  are  not  synonymous, 
they  are  often  used  interchangeably  in  the  literature  since  they  apply  to 
the  same  types  of  vehicles.     In  the  subsequent  discussion,   "continuous 
thrust"  is  used  to  describe  vehicles  which  thrust  for  a  major  portion  of 

the  transfer  and  "low  thrust"  refers  to  vehicles  which  have  acceleration 

_3 
levels  less  than  about  10      g 

6o 

This  study  is  concerned  primarily  with  those  propulsion  methods 
which  rely  upon  a  separate  energy  source*  for  the  generation  of  pro- 
pulsion energy.      With  few  exceptions,  these  are  low  thrust  devices 


>;<Underlined  words  are  defined  in  footnotes. 

separate  energy  source:     propulsion  system  characterized  by  a  power 
plant  which  is  independent  of  the  thrust-producing  mechanism. 


Where  it  is  necessary  to  be  more  specific,  an  on -board  nuclear  reactor 
is  assumed. 

We  shall  not  be  concerned  with  further  justifying  the  use  of  these 
devices,  nor  with  defining  their  regions  of  usefulness.     Such  questions 

9    *?    A    Pi 

are  well  covered  in  the  literature  .     The  goal  here  is  to  assume 

that  such  vehicles  will  exist  and  to  study  the  problem  of  guidance  ir- 
respective of  the  mission  or  of  the  particular  propulsion  method  used.. 

1.  1    The  Guidance  Problem 


The  presence  of  a  thrust  acceleration,  acting  over  a  significant 
portion  of  the  spacecraft  trajectory,  introduces  complexities  in  the 
analysis  of  spacecraft  motion  which  make  many  of  the  techniques  in 
common  use  for  ballistic  vehicles  inapplicable.     In  particular,  the  ele- 
gant conic  representations  can  be  used  only  when  the  thrust  may  be 
treated  as  a  small  disturbing  force.     In  general  this  is  the  case  only  in 
the  region  near  a  central  body.     Consequently  it  is  desirable  to  formu- 
late continuous-thrust  investigations  within  a  body  of  mathematics  hav- 
ing sufficient  applicability  to  treat  most  problems  of  interest.      The 
calculus  of  variations  and  the  concepts  of  optimal  control  theory  meet 
this  need. 

In  the  terminology  of  this  branch  of  applied  mathematics,  guidance 
of  continuous-thrust  spacecraft  is  the  problem  of  finding  a  control  pro- 
gram  which  transfers  the  spacecraft  between  two  given  points  in  space 
subject  to  constraints  on  the  velocity  at  both  terminal  points.     Since 
such  programs  are  not  unique,  it  is  possible  to  place  additional  con- 
straints on  the  solution  such  that  some  desired  quantity  (the  cost  func- 
tion) is  maximized  or  minimized.     Thus,  in  addition  to  satisfying  fixed 
boundary  conditions  at  both  launch  and  target  points,  solutions  may  be 
found  which  minimize  transfer  time,  or  propellant  consumed  or  some 


guidance:     as  used  in  this  thesis,  the  term  refers  to  the  process  of 
determining  a  control  program, 

control:     as  used  in  this  thesis,  the  term  refers  to  the  quantity  or  quan- 
tities (thrust  or  acceleration)  representing  the  propulsive  effort      It 
also  refers  to  the  mechanical  process  of  directing  the  propulsive  effort 
in  accordance  with  the  control  program. 


other  quantity  appropriate  to  the  mission,     This  thesis  is  concerned 

with  maximizing  only  the  final  mass.     At  the  present  time  no  method 

has  been  found  which  will  prove  rigorously  that  the  solutions  obtained 

in  inverse  square  (and  more  general)  gravitational  fields  satisfy  the 
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condition  for  an  absolute  maximum  of  the  final  mass         One  can  only 

show  that  the  solutions  produce  a  maximum  in  the  region  of  space  under 

consideration,  i   e     local  maxima,  and  then  demonstrate  from  physical 

reasoning  that  other  solutions  are  not  likely 

From  the  viewpoint  of  guidance,  initial  conditions  for  the  two -point 
boundary  value  problem  are  represented  by  the  present  state  of  the 
vehicle.     If  a  vehicle  can  be  flown  in  accordance  with  the  solution  that 
exists  at  the  instant  of  launch    only  one  solution  is  needed      Such  a  case 
is  a  problem  in  control,  not  in  guidance,     That  is,  the  desired  thrust 
program  is  known  (guidance)  and  the  magnitude  and  direction  of  thrust 
must  be  controlled  so  that  the   spacecraft    follows  the  desired  trajectory 
(control).     Even  with  "tight"  control  loops,,  however,  perturbations  in- 
evitably occur  which  cause  the  spacecraft  to  depart  from  tne  original 
optimized  trajectory      When  this  happens  a  new  control  program  must 
be  found  which  will  cause  transfer  from  the  present  state  to  the  final 
state  such  that  final  mass  is  maximized.     This  is  the  problem  we  desire 
to  solve.     It  follows,  that  if  a  solution  can  be  found  for  the  guidance  prob- 
lem, then  by  considering  the  launch  state  as  the  present  state,  the  op- 
timized trajectory  connecting  the  launch  point  and  target  point  may  be 
determined       This  latter  application  of  the  guidance  techniques  leads  to 
that  part  of  the  thesis  described  as  "trajectory  computation   ' 

In  the  preceding  paragraph  it  was  emphasized  that  the  present  state 
of  the  vehicle  serves  as  the  initial  condition  for  the  boundary  value  prob- 
lem.    It  is  necessary    therefore,  that  some  method  exist  for  determining 
the  state  of  the  vehicle  at   any  time.     In  this  thesis  the  unique  features 


state:     refers  to  the  seven  quantities  which  describe  the  vehicle  in  terms 
of  physical  variables 

present  state;  refers  to  the  vehicle  state  as  determined  by  the  space 
navigator,  i,  e  the  state  at  the  present  instant  of  time  In  this  thesis 
navigation  refers  to  the  process  of  determining  state  from  measure- 
ments 


of  a  continuous-thrust  vehicle  are  examined  from  the  viewpoint  of 
estimating  vehicle  state. 

Finally,  the  concept  of  constant  exhaust  power  is  examined  as  a 
criterion  for  low-thrust  transfers.     It  is  easily  shown  that  if  a  linear 
system  is  power -limited,  as  in  the  case  of  separately  powered  rockets, 
the  minimum  expenditure  of  energy  results  from  operation  at  maximum 
continuous  power    .      The  application  of  this  principle  has  been  used  ex- 
tensively throughout  the  literature  .     An  analogous  treatment  of 
thrust-limited  rockets  is  examined  here. 

1.  2    Prior  Studies 

Published  works  on  guidance  of  interplanetary  low-thrust  vehicles 
were  almost  nonexistent  prior  to  early  1963.     Miller       produced  a 
guidance  technique  for  cis-lunar  space  in  his  doctoral  research  in  1961 
This  technique  consists  of  spiralling  out  from  the  Earth  to  some  point 
from  which  the  vehicle  may  coast  to  the  vicinity  of  the  moon,  then 
matching  the  unique  velocity  vector,  corresponding  to  the  present  posi- 
tion of  the  vehicle,  which  will  result  in  achieving  the  target  point. 
Miller  showed  that  guidance  of  this  type  results  in  a  relatively  small 
fuel  penalty  for  lunar  missions. 

Friedlander       formulated  the  problem  in  the  classical  calculus  of 
variations  and  solved  the  adjoint  equations  in  two  dimensions  for  the 
sensitivity  coefficients  of  the  state  variables  along  an  optimized  tra- 
jectory to  Mars.     In  reference  12    he  suggests  a  linearized  solution 
which  minimizes  a  quadratic  function  of  the  control  variable  variation. 
The  solution  approaches  but  does  not  attain  the  control  program  for 
maximum  final  mass.     In  reference  13  Friedlander  applies  his  tech- 
niques to  a  vehicle  using  a  Snap- 8  power  source. 

14 
The  most  recent  work  is  that  of  Pfeiffer       who  applies  some  of  the 

newer  developments  in  the  theory  of  optimal  control  to  the  low-thrust 

guidance  problem.     Pfeiffer  solves  the  guidance  problem  and  minimizes 


power  limited:     refers  to  a  propulsion  system  characterized  primarily 
by  a  maximum  power  level. 

thrust  limited:     refers  to  a  propulsion  system  characterized  primarily 
by  a  maximum  thrust  level 


a  penalty  function  which  is  "equivalent"  to  a  quadratic  form  of  the  final 
state  error       The  control  program  produced  by  this  method  satisfies 
the  boundary  conditions  "as  closely  as  possible  using  the  penalty  func- 
tion  "    (the  quotation  marks  are  from  the  reference)    Pfeiffer's  method, 
however,  is  not  applicable  to  problems  where  certain  boundary  values 
are  fixed. 

Much  of  the  recent  work  in  optimal  control  theory  concerns  sys- 
tems which  have  mathematical  models  similar  to  those  for  low-thrust 
vehicles.      Many  of  the  ideas  developed  in  these  investigations  are 
directly  applicable  to  the  problem  considered  here.     Significant  con- 
tributions  are  attributable  to  Pontryagin      ,  Kalman       and  Breakwell 
who  have  formulated  existence  theorems  and  derived  necessary  and 
sufficient  conditions  for  optimal  trajectories       Breakwell  has  also 

contributed  important  work  in  specifying  the  form  of  solutions  with 
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constrained  control  vectors.     Athans,  Falb  and  LaCross       working 

together  and  individually  have  solved  many  special  cases  of  optimal 

trajectories  for  constrained  control  vectors. 

Fundamental  to  any  guidance  study  is  a  qualitative  knowledge  of  the 
trajectory  along  which  the  space  vehicle  is  to  travel      Quite  properly 
then,  the  earliest  work  in  low-thrust  propulsion  consisted  of  studies  of 
engine  characteristics  and  trajectory  characteristics.     Several  of  the 

earlier  studies  of  engine  characteristics  have  been  previously  refer- 

19 
enced.      To  those  must  be  added  the  contributions  of  Langmuir       and 

T      .      20 
Irving 

2  1 
In  the  area  of  trajectory  studies  Tsien       performed  some  of  the 

earliest  work  (1953).      This  was  followed  several  years  later  with  con- 
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tributions  by  Lawden      ,  Moeckel      ,  Melbourne      and  Zimmerman, 

25 
McKay  and  Rossa      .     The  problem  which  confronted  these  authors  was 

that  analytic  solutions  to  the  trajectory  problem  can  be  found  only  for 
linear  gravitational  fields.     For  central  force  fields  and  more  complex 
configurations,  solving  the  two-point  boundary  value  problem  was  a 
tedious  trial-and-error  procedure  requiring  the  use  of  high  speed  digit- 
al   computation      Satisfying  an  additional  constraint  for  an  optimized 
trajectory  was  tedious  and  time  consuming   even   with   high    speed 
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computers.     The  work  of  Bryson       and  others  in  the  late  1950's  and 
early  1960's  served  to  simplify  the  machine  procedures  so  that  com- 
putation of  optimal  trajectories  became  less  tedious  and  less  time 
consuming. 

1.  3    Thesis  Philosophy  and  the  Method  of  Approach 

Throughout  the  research  and  writing  of  this  thesis  the  author  has 
attempted  to  consider  the  low-thrust  guidance  problem  from  the  view- 
point of  a  space  navigator  who  is  responsible  for  the  safe  and  timely 
arrival  of  the  spacecraft,  at  the  target  point.     The  extrapolation  of  air- 
craft navigation  experience  into  space  navigation  is,  at  best,  a  hazard- 
ous undertaking;  however,  it  does  provide  a  basis  for  certain  decisions 
which  have  influenced  the  author's  approach  to  this  investigation.     The 
following  criteria  were  established  from  this  philosophy  and  have  been 
used  when  it  became  necessary  to  make  definite  assumptions. 

1)  The  mission  is  manned,  probably  utilizing  more  than  one 
vehicle,  each  of  which  is  manned  by  several  crewmen. 

2)  The  mission  duration  is  limited,  by  consideration  of  human 
tolerances. 

3)  The  spacecraft  configuration  and  the  mission  have  been  speci- 
fied.    Hopefully,  the  spacecraft  characteristics  are  optimum 
for  the  mission,  but  may  not  be. 

4)  Whatever  the  mission,  at  each  of  several  points  along  the  tra- 
jectory the  navigator  has  three  choices: 

a)  to  rendezvous  with  the  target  point  at  the  preplanned  time 
such  that  the  trajectory  minimizes  propellant  consumed. 

b)  to  rendezvous  with  the  target  point  utilizing  a  time  and 
trajectory  such  that  propellant  consumption  is  minimized. 

c)  to  rendezvous  with  the  target  point  in  minimum  time  utiliz- 
ing   the  available  propellant.     This  alternative  is  not 
treated  in  the  thesis. 

These  guidelines  establish  the  general  context  within  which  the  guidance 
problem  is  to  be  solved.     Let  us  now  proceed  to  examine  the  specific 
items  which  complicate  the  solution. 


1)  The  mathematical  theory  concerning  optimal  trajectories  dic- 
tates that  the  first  variation  of  the  optimized  quantity  must 
vanish  for  the  optimal  path  when  the  control  is  unconstrain- 
ed ;  As  a  consequence,  if  the  optimized  quantity  is  a 
state  variable    the  matrix  of  coefficients  relating  the  control 

variables  to  the  state  variables  is  singular  and  its  time  inte- 
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gral  along  the  optimal  trajectory  is  also  singular  .     The 

singularity  is  proven  in  Appendix  C.     Solutions  for  the  optimal 
control  usually  require  inversion  of  this  matrix.     The  prob- 
lem of  singular  matrices  is  handled  in  this  study  by  a  method 
of  deleting  certain  matrix  elements  which  create  the  singular- 
ity and  by  the  formation  of  a  new  matrix  which  can  be  inverted. 
The  deletion  method  is  an  important  part  of  the  thesis  and 
provides  a  general  method  for  treating  certain  singularities 
without  reformulating  the  problem 

2)  Optimal  trajectories  for  constrained  control  variables  often 
possess  discontinuities  in  the  first  derivatives  of  one  or  more 
state  variables  and  in  the  control  variables       This  mathemati- 
cal problem  is  handled  by  the  use  of  switching  functions  which 
are  continuous 

3)  Optimal  trajectories  for  constrained  control  variables  require 
periods  of  maximum  control  magnitude.     Therefore,  if  the 
maximum  propulsive  effort  is  required  for  the  optimum  tra- 
jectory, guidance  around  the  optimum  is  limited  to  changes  in 
thrust  direction  unless  reserve  propulsive  power  is  available 
for  guidance.     The  assumption  that  reserve  power  is  available 
is  used  for  this  study, 

It  was  stated  in  section  1    1  that  the  low-thrust  guidance  problem 
may  be  treated  as  a  two-point  boundary  value  problem  in  the  calculus 
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of  variations    "  Thus  techniques  of  the  calculus  of  variations 

constitute  a  primary  mathematical  tooL  One  of  these  techniques  is  the 
method  of  adjoints,  which  plays  a  fundamental  part  in  subsequent  chap- 
ters      The  author  has  borrowed  heavily  from  newer  theories  in  optimal 
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control      '        since  the  state  space  formulations   widely  used   in   the 


literature  of  that  field  are  applicable  to  low -thrust  guidance.     One  of 
the  more  useful  tools  in  optimal  control  theory  is  associated  with 

1    c  1   ft 

Pontryagin       although  other  authors  have  used  the  same  principle 

A  derivation  of  Pontryagin' s  maximum  principle  is  outlined  in 
Appendix  D  for  convenience  of  the  reader. 

To  facilitate  notation  and  to  preclude  the  possibility  of  fundamental 
notions  becoming  obscured  by  the  quantity  of  algebraic  detail,  matrix 
notation  and  the  ideas  of  matrix  calculus  are  used  throughout..     Finally, 
to  test  the  thesis,  numerical  analysis  and  an  IBM  7094  were  employed. 

1.  4   Relationship  to  Prior  Studies 

In  the  research  preceding  this  study,  the  author  began  with  the 

1  D  R   9  4- 

formulations  of  Friedlander       and  Melbourne  ,  and  attempted  to 

extend  their  ideas  into  areas  of  more  general  application  and  to  find 
solutions  to  the  guidance  problem  which  were  useful  from  the  space 
navigator's  viewpoint.     One  of  the  fundamental  considerations  was  to 
find  control  laws  which  optimize  propellant  consumption.     Cost  func- 
tions which  produce  near-optimal  propellant  consumption  were  rejected. 

The  idealized  formulations  for  the  separately  powered  rocket  re- 
quire a  wide  range  of  thrust  and  of  specific  impulse  as  the  vehicle 
traverses  its  trajectory.     Current  technology  indicates  that  variable - 
specific -impulse  thrusters  will  not  be  available  in  the  forseeable 
future,  at  least  for  electrostatic  vehicles.     To  satisfy  this  engineering 
restriction,  investigators  have  continued  to  use  the  power-limited  for- 
mulations but  approximate  the  optimal  thrust  magnitude  programs  with 
regions  of  constant  specific  impulse  where  relatively  low  values  of 
specific  impulse  are  required  and  with  coast  elsewhere.     It  is  shown 
in  this  thesis  that  the  so-called  "bang-bang"  control  used  to  satisfy  the 
engineering  restriction  can  be  derived  by  abandoning  the  concept  of 
power-limited  thrusting. 

In  section  1.  1  the  necessity  of  estimating  the  state  is  discussed. 
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The  works  of  Battin      ,  Stern      ,  and  Potter  and  Stern       serve  as  the 


electrostatic  (propulsion)  :     A  propulsion  method  depending  upon  the 
acceleration  of  charged  particles  through  an  electrostatic  field. 


starting  point  for  extending,,  to  the  low-thrust  case,  the  navigation  tech- 
niques (estimate  of  state)  which  have  been  developed  for  ballistic  vehicles. 

Finally,  using  the  guidance  solutions  an  iterative  technique  for  com- 
puting low-thrust  trajectories  is  developed.     Recent  investigations  in  the 
field  of  trajectory  computation  appear  to  rely  heavily  on  steepest  ascent 
techniques       These  techniques  require  rather  complex  programs,  care- 
ful handling  and  suffer  from  the  problem  of  slow  convergence  as  the 
optimal  trajectory  is  approached.     The  technique  suggested  in  this 
study  appears  to  be  faster  even  for  discontinuous  control  variables, 

1,  5   Summary  of  Thesis  Objectives 

In  closing  this  introductory  chapter  it  seems  appropriate  to  provide 
a  sketch  of  the  objectives  of  subsequent  discussions       Briefly,  the  ob- 
jectives of  this  thesis  are:     1)  to  present  a  linear,  noniterative  method 
for  deriving  a  control  program  for  guidance  of  low-thrust  space  vehicles 
with  respect  to  a  propellant- optimal  reference;  2)  to  show  that  the  com- 
puted control  programs  satisfy  the  necessary  conditions  for  an  optimum 
3)  to  show  that  the  method  may  be  extended  to  trajectory  computation  by 
successive  iteration  of  the  guidance  solution,  4)  to  examine  a  method  of 
estimating  the  state  of  continuous- thrust  vehicles  and  5)  to  examine  the 
concepts  of  power- limited  and  thrust-limited  vehicles  from  the  view- 
point of  guidance. 

The  fundamental  argument  of  the  thesis  may  be  extracted  from 
these  several  objectives       It  is:      'There  exists  a  linear  method  which 
produces  a  propellant-optimal  control  program  in  a  noniterative  form 
for  guidance  of  power-limited  and  thrust-limited  space  vehicles,  and 
which  provides  a  simple,  rapidly  converging  iterative  technique  for  com- 
puting propellant-optimal  trajectories'  ' 


propellant-optimal:     refers  to  a  trajectory  which  results  in  minimum 
propellant  usage. 


CHAPTER  II 
THE  LOW-THRUST  SPACE  VEHICLE 

2.  1   Summary  of  Chapter  II 

In  this  chapter  the  parameters  which  characterize  low-thrust  vehicle 
performance  are  derived  and  discussed.     The  constant -exhaust-power 
concept,  which  has  been  widely  used  in  the  past,  is  shown  to  be  the  op- 
timum method  of  engine  control.     The  propellant  cost  accruing  from 
engineering  restrictions  on  variable  specific  impulse  is  computed  for 
field-free  space  by  analytic  methods.     Both  methods  of  engine  control, 
the  ideal  and  the  practical,  are  discussed  from  the  viewpoint  of  guidance, 
Finally,  the  results  of  a  numerical  study,  which  confirm  the  derivations, 
are  presented. 

2.  2    The  Use  of  a  Separate  Energy  Source 

One  may  verify  from  momentum  considerations  that  the  instantane- 
ous force  exerted  on  a  space  vehicle  by  its  exhaust  stream  is: 

f_=  mc  (2-1) 

where  f_is  the  force  vector,  c  is  the  oppositely  directed  exhaust  velocity 
and  m  is  the  mass  rate  of  the  vehicle,      (-m  is  the  propellant  flow  rate    ) 

Since,  the  force  magnitude  may  be  held  constant  if  m  and  c  are 
varied  inversely,  the  selection  of  a  high  exhaust  velocity  and  low  pro- 
pellant flow  tends  to  reduce  the  total  amount  of  propellant  required  for 
a  given  impulse.     Unfortunately,  processes  which  use  the  products  of 
combustion  as  the  propellant  are  unable  to  produce,  simultaneously,  the 
low  flow  rates  and  high  exhaust  velocities  required  for  many  interesting 
missions.     If  we  use  specific  impulse,  I      ,  as  a  measure  of  the  engine 
performance,  where 

c                 f 
I       =  ___    =  1 (2-2) 

SP        «.  ,     '  \ 

g0  (-m)  gQ 
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and  g     is  the  acceleration  due  to  gravity  at  the  Earth's  surface,  then 

chemical  systems  are  limited  to  values  of  specific  impulse  under  600 

5 
seconds  and  direct  nuclear  systems  to  values  of  about  1600  seconds 

However,  if  a  separate  energy  source  is  available,  the  source  energy 

may  be  converted  into  electrical  energy  and  used  to  accelerate  charged 

propellant  particles.      Values  of  specific  impulse  in  the  range  of  3000 

to  20,000  seconds  appear  attainable  in  this  way 

Several  energy  sources  and  conversion  processes  are  under  inves- 
tigation.    These  may  be  divided  into  the  two  broad  categories  of  direct 
or  indirect  energy  conversion..     Direct  methods  may  be  characterized 
by  the  absence  of  a  mechanical  phase  in  the  conversion  process.     Power 
from  solar  cells  is  cne  example  of  this  type.     Direct  production  of  elec- 
trical energy  in  a  nuclear  reactor  is  another.     A  proposal  for  this  latter 
method  is  discussed  in  reference  33    however  it    has  not  been  proven  in 
the  laboratory. 

Indirect  conversion  of  the  source  energy  into  electrical  energy 
appears,  at  the  present  time,  to  be  more  realistic  for  large  manned 
spacecraft  which  have  power  requirements  in  the  range  of  several 
megawatts       The  process  generally  considered  most  promising  uses  a 
nuclear  source  to  power  turbomachinery  for  the  production  of  electrical 
energy.     A  block  diagram  of  this  type  system  is  presented  in  Figure  2-a. 

In  this  report,  an  on-board  energy  source  is  explicitly  assumed 
although  much  of  the  guidance  analysis  is  independent  of  such  considera- 
tions.     This  assumption  permits  the  power  availability  to  be  independent 
of  the  trajectory.     This  is  not  the  case,  for  example,  when  solar  energy 
collectors  are  used  since  for  a  given  collector  area  the  power  availa- 
bility varies  inversely  with  the  square  of  the  distance  to  the  sun. 

2,  3    The  Constant-Power  Concept 

Separate  energy  sources  are  most  often  described  as  power -limited 
devices       That  is,  their  rate  of  energy  conversion  is  the  design  criterion. 
If  the  energy  produced  by  the  source  is  all  converted  to  propulsive  en- 
ergy then  the  power  in  the  exhaust  stream  is  given  by 
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Fig.  2-a.    Schematic  diagram  of  separately  powered  rocket. 


P 


-m  c 


(2-3) 


(The  efficiency  of  the  conversion  process  will  be  introduced  in  a  sub- 
sequent discussion.  ) 

It  is  well  known  that  for  a  large  class  of  problems  in  linear  fields, 
minimum  total  energy  expenditure  results  from  operation  at  maximum 
continuous  power.      This  result  may  be  obtained  for  conservative  fields 
in  general  by  use  of  variational  techniques  as  shown  in  reference  7.    The 
usual  assumption  in  low-thrust  investigations,  that  exhaust  power  is  a 
constant  and  equal  to  its  maximum  value,  is  therefore,  quite  valid  in  the 
idealized  problem.      This  result  is  used  to  establish  parameters  for  the 
separately  powered  rocket. 
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Combining  equation  (2-3)  with  the  relations 

f  =  -mc  (2-4) 

f  =  ma  (2-5) 

one  obtains 


2 

■m  a 


2  o 

m  2p 


(2-6) 


where  a  is  the  thrust  acceleration  of  the  vehicle.     Integration  of  equa 
tion  (2-6)  over  the  flight  path  yields 

i  i  l       r     2 


f    a     dt  (2-7) 


mr  m  2p      0 

f  o  ^ 

where  the  subscripts  indicate  initial  and  final  times      Equation  (2-7)  is 
more  conveniently  expressed  in  terms  of  the  mass  ratio. 

tf 
m  r        o 

MR  -   1    =    —      J      a     dt  (2-8) 

2p      0 

Three  parameters  will  now  be  defined 

J  =    i     /    a2  dt  (2-9) 

Z    0 

m 
a    =   —EL.  (2-10) 

P 

(3     =    — —  (2-11) 

m 
o 

J   is  the  well  known  acceleration  integral,  a  is  the  specific  mass  of  the 
propulsion  system  and  must  include  the  total  efficiency  of  energy  con- 
version when  p  is  the  desired  exhaust  power,  (3  is  the  propulsive  system 
mass  fraction  and  m^  is  the  propulsive  system  mass.     Substituting 
these  parameters  into  equation  (2-8)  one  obtains 

MR    =     1  +  -    J  (2-12) 
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For  a  spacecraft  with  fixed  parameters  and  a  specified  mission, 
the  minimum  mass  ratio  and  thus  the  minimum  propellant  usage  is 
achieved  for  trajectories  which  make  J  a  minimum. 

The  objective  of  trajectory  computation  and  of  guidance  is  to  find 
a  path  which  accomplishes  the  mission  and  minimizes  J. 

2.  4    Limitation  of  the  Constant-Power  Concept 

Using  the  equations  of  section  2.  3,  the  thrust  acceleration  may  be 
expressed  as 


_    2p 


cm 


(2-13) 


Figure  2-b  shows  a  plot  of  the  thrust  acceleration  magnitude  during  a 
fast  transfer  to  Mars  for  which  J  is  a  minimum.      The  absolute  value 
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Fig.  2-b.    Optimum  acceleration  level  for  variable-specific-impulse  Earth-Mars  transfer. 
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of  acceleration  has  a  variation  over  two  orders  of  magnitude,      From 
equation  (2-13)  it  is  apparent  that  the  exhaust  velocity,  (i„  e.  ,  the  specif- 
ic   impulse)  must  have  a  similar  range  of  values  since  m  is  monotoni- 
cally  decreasing      However    the  current  thinking  of  investigators  who 
are  studying  propulsive  devices  indicates  that  in  a  given  thrust  pro- 
ducer, only  a  very  small  range  of  exhaust  velocity  variation  will  be 
feasible  unless  major  advances  in  technology  are  forthcoming.     The 
thrust  producer  can  be  designed  for  one  particular  exhaust  velocity. 
Operation  at  other  than  the  design  point  tends  to  decrease  efficiency 
and  shorten,  markedly,  the  useful  life  of  the  device. 

19 
This  restriction  on  exhaust  velocity  is  treated  in  the  literature 

as  a  departure  from  the  ideal  conditions  of  constant  exhaust  power. 
The  author  treats  this  restriction  from  a  slightly  different  viewpoint. 
In  the  succeeding  section  parametric  equations  for  a  constant-specific- 
impulse  rocket  will  be  derived  which  are  analogous  to  equations  (2-8) 
through  (2-12)       With  these  parametric  equations  one  may  compare 
directly  the  performance  of  a  rocket  with  specified  maximum  power, 
under  the  two  modes  of  propulsive  control,  L  e.   constant-specific-im- 
pulse (CSI)  or  variable-specific-impulse  (VS1) 

2,  5    The  Constant-Specific-Impulse  Concept 


The  propellant- optimal  thrust  program  for  a  CSI  rocket  is  shown 
subsequently  to  have  two  magnitudes,  maximum  thrust  or  zero  thrust. 
In  order  to  allow  for  maneuvering  and  for  departures  from  the  refer- 
ence trajectory,  reserve  thrust  capability  is  required.     (For  VSI  rock- 
ets this  problem  does  not  occur  since  thrust  is  continuously  variable.  ) 
A  method  for  controlling  thrust  magnitude  which  appears  quite  attrac- 
tive and  is  compatable  with  thruster  design  is  to  build  up  the  propulsion 
unit  from  a  large  number  of  small  thruster  nozzles  of  constant  specific 

impulse.     Proposals  of  this  type  appear  in  the  literature  but  are  not 

35 
analyzed  in  detail  Control  of  thrust  magnitude  is  obtained  by  con- 

trolling the  number  of  thrusters  in  operation 

If  a  nuclear  energy  source  is  used  the  power  range  is  not  continu- 
ously variable  from  zero  to  maximum  power  and  may  be  limited  in  the 
number  of  times  that  stopping  and  starting  is  permitted      However,  the 
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energy  source  and  the  thruster  units  can  be  operated  at  their  respective 
optimums  by  placing  an  external  high  temperature  resistor  in  the  sys- 
tem which  will  radiate  energy  in  excess  of  that  required  by  the  thrusters. 
Since  the  nuclear  fuel  constitutes  only  a  very  small  fraction  of  vehicle 
weight,  this  procedure  has  negligible  effect  on  the  final  mass  of  the 
vehicle,      Figure  2 -a  shows  the  propulsive  system  schematic  diagram 

For  a  propulsion  unit  constructed  in  this  manner,  the  total  thrust 
of  the  vehicle  is  some  fraction  of  the  total  thrust  available       That  is 

(2-14) 

(2-14a) 

(2-14b) 


(2-15) 


f 

- 

-nm   c 
o 

f 

o 

= 

-m   c 
o 

m 

= 

nm 
o 

nf 
o 

= 

-nm   c 
o 

m 

m 

P  =  npQ 

= 

2 
-nm   c 
o 

(2-16) 


where  n  is  the  fraction  of  thrusters  in  use  and  f  ,   -m     and  p     are  the 

o  o  ^o 

maximum  values  of  thrust,  propellant  flow  rate  and  exhaust  power., 
respectively. 

By  eliminating  the  exhaust  velocity,  equations  (2-14)  through  (2-16) 
may  be  manipulated  into  the  form 

-nm  f    a 

=  — - (2-17) 


2  o 

m  2pm 

Integrating  and  expressing  the  result  in  terms  of  the  mass  ratio,  one 
obtains 


m 


tf    f 


MR  -   1  =  — —     /      —     a  dt  (2-18) 


2p        o  m 

^o 


Observe  that  the  factor  f  /m  is  the  maximum  acceleration  possible 

o' 

for  a  given  mass.     Therefore  the  acceleration  integral  for  a  thrust- 
limited  vehicle  is  proportional  to  the  integral  over  the  thrusting  time  of 
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the  product  of  the  maximum  instantaneous  acceleration  and  the  instan- 
taneous acceleration.     Designate  this  integral  J*  and  observe  that  for 
a  vehicle  of  specified  power  and  initial  mass,  the  mass  ratio  is  a  min- 
imum when  J*  is  a  minimum 

J*     =    -      f    a  a  dt  (2-19) 

2    0J        max 

MR  -   1  =    -     J*  (2-20) 

P 

Equations  (2-19)  and  (2-20)  correspond     to  equations  (2-9)  and 
(2-12)  respectively,  and  permit  a  direct  comparison  of  the  performance 
of  a  hypothetical  vehicle  with  given  power  and  given  initial  mass  when 
controlled  as  a  VSI  vehicle  and  then  as  a  CSI  vehicle. 

The  comparison  is  made  by  finding  an  optimal  trajectory  for  a 
specified  mission  for  each  type  of  control  and  then  comparing  the  mass 
ratios  or  the  acceleration  integrals  for  the  two  cases.     This  has  been 
done  for  a  simulated  mission  to  Mars  using  techniques  described  in 
later  chapters,  and  for  the  simple  case  of  a  transfer  in  field-free  space, 

2,  6    Comparison  of  CSI  and  VSI  Control 


From  a  qualitative  point  of  view,  finding  an  optimal  trajectory  for 
the  VSI  machine  consists  of  finding  that  trajectory  which  minimizes 
propellant  expenditure  for  a  given  constant  exhaust  power  but  with  un- 
constrained exhaust  velocity.     In  the  case  of  the  CSI  vehicle,  the  prob- 
lem is  that  of  finding  the  trajectory  which  minimizes  the  propellant 
expenditure  subject  to  a  given  maximum  exhaust  power  and  a  given 
constant  exhaust  velocity 

If  both  vehicles  have  the  same  maximum  exhaust  power  and  both 
control  programs  optimize  propellant  consumption  then  constant-specific- 
impulse    is    but    a    limiting  case  of  variable-specific-impulse  and 
cannot  result  in  less  propellant  usage.     Therefore  the  acceleration 
integral  for  a  VSI  transfer  can  be  used  as  a  reference  for  assessing 
the  additional  propellant  cost  of  CSI  transfers 
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In  order  to  compare  the  two  vehicles  for  a  transfer  in  the  gravita- 
tional field  it  is  necessary  to  use  machine  computation.     However, they 
may  be  compared  analytically  in  field-free  space.     The  results  are 
shown  to  provide  reasonable  approximations  for  the  values  obtained  in 
the  gravitational  field. 

Assume  that  in  field-free  space  it  is  desired  to  traverse  the  dis- 
tance L,  in  the  time  T  such  that  the  rocket  begins  and  ends  at  rest.     This 

34 
problem  is  solved  in  the  literature  for  VSI  spacecraft      .     The  solution 

is  rederived  in  Appendix  A  and  the  solution  for  the  CSI  spacecraft  is 

also  derived. 


The  results  for  variable-specific-impulse  thrusting  are: 

2 
J 


6L' 


min 


T 

6L 

lo    "  T2 


3 


(2-21) 


(2-22) 


where  a     is  the  initial  acceleration.     The  optimum  acceleration  pro- 
gram starts  at  a     and  decreases  linearly  with  time  such  that  a(T)  =  -  a   . 

&>  0  j  \     i  0 

(See  Figure  2-c.) 


Fig.  2-c.    Optimum  acceleration  for  constant  power  rocket  in  field-free  space. 
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For  a  normalized  mass  (i.  e     m     is  unity),,  equation  (2-12)  may  be 

o 


written  in  this  case  as: 


Define  the  parameter  R. 


2 
MR  =  ^-    +    1  (2-23) 

PTJ 


R  =  —  (2-24) 

P 

MR  =  —    +    1  (2-25) 

3 

Therefore  in  this  case 

2 
R  =  ML.  (2-26) 

Equation  (2-25)  is  plotted  in  Figure  2-d.     R  is  an  excellent  measure 
of  the  difficulty  of  a  mission  in  field-free  space.     For  example:     A 
large  transfer  distance,  a  short  transfer  time  and  a  small  energy 
source  would  produce  a  large  value  of  R      Such  conditions  indeed  repre- 
sent a  difficult  transfer,, 

It  should  be  noted  that  the  mass  ratio  is  linear  in  R,      Thus  VSI 
transfers  characterized  by  a  finite  R  require  finite  mass  ratios, 

For  a  constant-specific-impulse  vehicle  of  identical  power  to  per- 
form this  mission  requires  a  specification  of  exhaust  velocity  or,  more 
conveniently,  initial  acceleration.     In  Figure  2-d,  the  mass  ratio  is 
plotted  as  a  function  of  R  for  two  values  of  initial  acceleration.      The 
first  corresponds  to  a     =  — ~-    and  is  asymptotic  to  R  -  20,     The  second, 

which  corresponds  to  an  optimum  value  of  initial  acceleration  for  CSI 
machines  is  exponential  in  R 

The  optimum  initial  acceleration  for  CSI  vehicles  is 

aQ  =    2     (in  MR)  aQ  (2-27) 

(opt)CSI        R  (OP^VSI 

For  small  values  of  R  such  that  the  mass  ratio  is  near  unity,  equation 
(2-27)  may  be  simplified  to 
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a  =  —    a  (2-28) 

°(opt)csl        4        °(opt)VSI 

The  transcendental  equations  for  this  example  are  derived  and  solved 
in  Appendix  A. 

The  conclusions  that  may  be  drawn  from  the  field-free  space 
example  are: 

1)  Control  of  thrust  by  varying  specific  impulse  but  with  con- 
stant exhaust  power  is  the  propellant-optimal  method  of  con- 
trol. 

2)  For  missions  represented  by  small  values  of  the  parameter 
R,  (less  than  5)  the  propellant  penalty  for  using  an  optimized 
constant  specific  impulse  is  less  than  15%, 

3)  For  large  values  of  R  the  propellant  cost  of  constant-specific- 
impulse  is  very  large,  even  for  an  optimized  CSI 

4)  For  sufficiently  long  transfer  times  the  value  of  R  may  be 
decreased  so  that  the  propellant  penalty  of  constant-specific- 
impulse  is  very  smalL      (As  will  be  subsequently  shown,  an 
optimized  one  way  trip  to  Mars  corresponds  to  an  R  of  approx- 
imately 1). 

It  is  not  readily  apparent  that  the  results  derived  for  field-free 

space  are  directly  applicable  to  transfers  in  the  gravitational  fields  of 

34 
the  solar  system.     However,  Melbourne  and  Sauer       have  computed 

the  VSI  acceleration  requirements  for  a  number  of  interplanetary  trans- 
fers and  found  the  values  to  be  in  surprisingly  good  agreement  with  the 
field-free  space  analysis.     Therefore,  it  seems  valid  to  argue  on  the 
basis  of  field-free  space  when  comparing  CSI  and  VSI  systems.     The 
numerical  results  in  this  study  support  this  conclusion. 

The  parameter  R,  which  includes  the  effects  of  both  mission  re- 
quirements and  spacecraft  power,  was  found  to  be  more  useful  in  the 
comparison  of  CSI  and  VSI  systems  than  the  acceleration  integral  alone. 
Some  significant  facts  which  are  evident  when  R  is  used  as  a  parameter 
but  which  are  missed  when  the  acceleration  integral  alone  is  used  are: 
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1)  increasing  the  spacecraft  power  for  a  particular  mission 
reduces  the  propellant  penalty  of  CSI  operation, 

2)  the  optimum  initial  acceleration,  and  thus  the  optimum  ex- 
haust velocity  for  CSI  vehicles  may  be  determined  to  reason- 
able accuracy  from  analysis  of  the  VSI  spacecraft  using  equa- 
tion (2-27). 

2.  7    Maximizing  the  Payload 

In  section  2.  6,  the  problem  of  maximizing  the  mass  ratio  is  con- 
sidered. In  this  section  the  energy  source  which  maximizes  the  pay- 
load  for  a  given  mission  and  mode  of  operation  is  computed. 

Both  CSI  and  VSI  vehicles  are  characterized  by  the  general  equa- 
tion: 

MR  =    -     J  +  1  (2-29) 

P 

m 
where  the  term  ~-  J  may  be  written  alternatively  as  J. 


If  the  initial  mass  of  the  vehicle  is  considered  to  consist  of  only 
three  parts:    payload,  propellant  and  power  source,  then  the  mass  dis- 
tribution may  be  written  as 

(2-30) 


(2-31) 


where  

m 
o 


mL 

=     1  -  (3 

-  (1  - 

-      '     ) 

m 
o 

MR 

mL 
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(3 

m 
o 

MR 

is  the  payload  fr 

action. 

Substituting  equation  (2-29)  into  equation  (2-31),  an  expression  for 
payload  fraction  is  obtained  in  terms  of  a,  (3,  and  J. 


-    1  (2-32) 
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Differentiating  equation  (2-32)  with  respect  to  (3  and  setting  the  deriva- 
tive equal  to  zero  results  in  the  optimal  mass  distribution  for  either 
CSI  or  VSI  space  vehicles.     The  result  is: 

—  J  =  (1  -    V~^J)2  (2-33) 


m 
o 


opt 


V2 


Popt  =  (a  J)  '       "  aJ  (2-34) 


(2-35) 


m 
m 


here  — —  is  the  propellant  fraction. 


o 

From  the  alternative  forms  of  equation  (2-29)  and  equation  (2-34) 
the  optimum  exhaust  power  for  a  space  mission  is  obtained. 

Popt     =    Po    =     [J\     l}2     _  j  (2_36) 


a  m  v  a 

o 

For  the  numerical  work  in  this  thesis  an  optimistic,  but  not  un- 
reasonable, value  of  a  was  selected.     In  all  subsequent  computations, 
an  a  of  10  kilograms  per  kilowatt  is  assumed.     From  the  digital  com- 
puter studies  of  a  150  day  Earth-Mars  one  way  transfer,  a  value  of  J 
was  obtained  for  the  sizing  relationships.     This  value  of  J  was  multi- 
plied by  the  factor  4  and  rounded  off  to  provide  a  more  realistic  value 
for  a  round  trip  to  Mars.     The  use  of  an  approximation  is  justified  in 
that  the  purpose  here  is  to  provide  a  reasonable  value  of  power  for  com- 
paring VSI  and  CSI  operations  without  placing  undue  emphasis  on  opti- 
mizing the  round  trip  mass  distribution       That  is,  the  author  is  more 
concerned  with  optimal  guidance  of  a  space  vehicle  of  given  power  and 
mode  of  operation  than  in  determining  whether  the  spacecraft  is  the 
best  vehicle  for  the  mission       Further,  the  preceding  method  of  maxi- 
mizing payload,  although  widely  used  in  low-thrust  studies,  does  not 
adequately  treat  a  round  trip  mission  comprised  of  several  phases  such 
as  escape  from  Earth,  transfer,  capture  at  the  target,  a  waiting  period, 
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and  similar  return  phases.  Each  of  the  phases  may  impose  particular 
requirements  which  are  not  compatible  with  the  other  phases.  For 
example,  it  may  be  impractical  to  use  100  or  more  days  to  effect  a 
spiral  escape  even  though  the  optimized  solution  requires  it.  Conse- 
quently the  derivation  in  this  section  must  be  considered  as  only  an 
approximate  method  of  maximizing  payload  for  realistic  round  trip 
mission  planning. 

—  fi 
The  value  of  J  selected  was  10      .     The  units  are  (astronomical 

2  3 

units)     per  (day)    .     Inserting  these  values  of  a  and  J  (in  compatable 

units)  into  the  equation  (2-36),  a  value  for  p   /m     was  computed  for  use 
in  the  numerical  work. 

— -      =      0.  0242  kw/kg 

m  (2-37) 

o  v  ' 

=      0.  6988  X  10"6  (A.  u.  )2/(day)3 

Using  the  value  obtained  in  equation  (2-37),  values  of  MR  were 
extracted  from  computer  results  and  compared  with  the  field-free 
space  values  for  both  CSI  and  VSI  engine  control.     The  results  of  this 
comparison  are  plotted  on  the  curves  of  Figure  2-d.     Agreement  is 
good  for  the  cases  tested  and  warrants  further  investigation  for  dif- 
ferent values  of  the  parameter  R. 

Values  of  initial  acceleration  and  specific  impulse  which  are  char- 
acteristic of  the  values  obtained  in  the  numerical  studies  are 

a   la     =  1.  2  X  10"4  (2-38) 

o/&o 

I        =  4000  seconds  (2-39) 

sp  ' 
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CHAPTER  III 
LINEAR  GUIDANCE 

3.  1   Summary  of  Chapter  III 

The  guidance  problem  for  continuous -thrust  vehicles  is  formulated 
in  this  chapter  using  perturbation  techniques  and  the  theory  of  optimum 
control.     A  reference  trajectory  which  satisfies  mission  requirements 
is  assumed  to  exist  and  its  associated  control  program  is  known      The 
consequences  of  allowing  the  reference  trajectory  to  be  an  optimal  tra- 
jectory are  examined.      Linearized  equations  of  state  are  solved  to- 
gether with  their  adjoint  set  to  produce  a  state  transition  equation 
applicable  for  small  perturbations  around  the  reference  trajectory 
The  state  transition  equation  is  then  shown  to  be  suitable  as  a  guidance 
equation.     Solutions  for  the  guidance  equation  are  derived  for  fixed  - 
time-of-arrival  (FTA)  guidance.      Variable-time-of-arrival  (VTA)  guid- 
ance is  solved  for  certain  restricted  cases. 

3.  2    General  Remarks 


In  terms  of  the  discussion  in  Chapter  II,  the  ensuing  derivation  of 
linear  guidance  may  be  characterized  as  a  method  of  minimizing  the 
acceleration  integral,  J,  between  the  present  position  and  the  target 
position.      For  this  purpose  it  is  convenient  to  abandon  temporarily  the 
engineering  aspects  of  low-thrust  transfers  and  to  consider  the  prob- 
lem in  terms  of  the  calculus  of  variations  and  the  theory  of  optimum 
control. 

The  guidance  schemes  suggested  for  interplanetary  vehicles  fol- 
lowing ballistic  trajectories  usually  require  that  midcourse  corrections 


state  transition  equation;.  An  equation  which  expresses  the  state  of  the 
vehicle  at  any  given  time  in  terms  of  1)  the  state  at  any  other  time  and 
2)  the  control  existing  between  the  two  times 

guidance  equation:     An  equation,  possessing  a  solution  which  satisfies 
the  mission  requirements  for  the  spacecraft.      The  solution  of  a  guid- 
ance equation  is  often  termed  the  "control  law"  or  "control  program" 
for  the  physical  guidance  system. 


25 


minimize,  subject  to  one  or  more  constraints,  the  position  and  veloc- 
ity   variations  at  the  destination  point. 

The  characteristics  of  a  continuous-thrust  vehicle  permit  a  guid- 
ance scheme  which  will  null  the  position  and  velocity  variations  at  the 
destination,  provided  guidance  is  begun  sufficiently  early  in  the  trans- 
fer.    For  all  practical  vehicles,  if  guidance  is  not  initiated  early  in  the 
transfer  the  low-thrust  vehicle  will  arrive  at  some  point,  beyond  which 
insufficient  power  is  available  to  complete  the  mission.     The  vehicle 
is  then  said  to  be  "in  extremis".     This  boundary  may  be  visualized  as 
a  conical  surface  surrounding  the  reference  trajectory  with  its  apex 
at  the  target  point.      The  "distance"  in  phase  space  from  the  trajectory 
to  the  surface  depends  upon  the  excess  power  allowed  for  guidance. 
The  problem  of  determining  the  optimum  power  and  propellant  allow- 
ance for  guidance  is  not  studied  in  this  report;  however,  the  techniques 
of  this  chapter  are  readily  applicable  to  such  investigations  and  are 
recommended  as  a  tool  for  future  work.     In  this  chapter  it  is  assumed 
that  the  reserve  power  is  adequate  for  satisfying  the  solution  to  the 
guidance  equation. 

3.  3    Formulation  of  the  Problem 

The  basis  for  the  subsequent  analysis  is  that  a  vehicle  is  in  transit 
betwe?n  two  planets.      The  control  program  in  use  corresponds  to  some 
known  reference  trajectory.     Measurement  of  vehicle  state  indicates 
that  the  vehicle  is  off  the  reference  trajectory  by  some  small  amount 

The  first  problem  of  interest  asks  the  question:     "What  is  the 
effect  of  the  state  variation  at  the  time  of  the  measurement  on  the  ve- 
hicle state  at  any  future  time?"    The  next  obvious  question  is:    "What 
is  the  new  control  program  which  will  cause  the  vehicle  to  satisfy 
mission  requirements?"    The  third  is:     "If  there  is  more  than  one  con- 
trol program  available,  what  criterion  should  be  used  for  selecting  one 
of  them?" 

The  third  question  may  be  disposed  of  immediately.     The  thesis 
is  concerned  with  propellant-optimal  guidance.     Therefore  the  choice 
between  controls  is  on  the  basis  of  minimum  propellant  consumption. 
From  Chapter  II,  a  criterion  which  satisfies  this  objective  is  the 
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acceleration  integral.  There  may  be  other  criteria  which  are  equally 
as  good,  however  J  was  chosen  because  it  is  easy  to  manipulate  and  is 
easily  interpreted  in  the  physical  problem 

The  second  question  requires  that  the  mission  be  defined.      For 
this  purpose  it  is  sufficient  to  require  that  the  final  position  and  veloc- 
ity   of  the  vehicle  have  certain  well  defined  values. 

Depending  upon  the  particular  mission,  the  above  guidance  object- 
ives may  be  satisfied  at  the  preplanned  arrival  time  (FTA  guidance) 
or  at  some  time  different  from  the  planned  time  but  consistent  with  the 
mission  (VTA  guidance).     The  FTA  problem  is  readily  solved  by  the 
linear  methods  of  this  chapter.     The  VTA  propellant-optimal  problem 
is  considerably  more  difficult      A  method  of  obtaining  restricted  solu- 
tions by  linear  methods  is  presented  and  discussed. 

3.  4   Selection  of  State  Variables 


The  differential  equations  of  motion  for  celestial  bodies  admit  six 
constants  of  integration.      The  selection  of  the  six  quantities  to  repre- 
sent this  motion  is  to  some  extent  arbitrary.     To  be  consistent  with 
the  definition  of  mission  requirements  in  this  investigation,  and  be- 
cause they  are  convenient,  the  three  components  of  position  and  the 
three  components  of  velocity  are  chosen      The  convenience  is  due  to 
their  facility  for  visualization  and  their  relative  facility  for  physical 
measurement, 

Since  the  conservation  of  propellant  is  important,  an  additional 
differential  equation  describing  the  mass  change  due  to  propellant  flow 
must  be  included  for  certain  cases.      Thus  seven  independent,  but 
coupled,  variables  are  sufficient  to  describe  the  state  of  the  vehicle 
at  any  time,      The  state  will  be  subsequently  written  as  a  vector  which 
has  components  of  position,  velocity  and  mass. 


(3-1) 


The  second  order  differential  equations  of  motion  will  be  rewritten  as 
first  order  equations  to  conform  to  the  above  definition  of  state. 
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For  this  thesis,  only  central  force  fields  are  considered      It  is 
subsequently  shown  that  extension  to  more  general  gravitational  fields 
does  not  invalidate  the  theory  but  does  introduce  computational  com- 
plexities which  serve  only  to  obscure  the  physical  concept  if  introduced 
prematurely. 

With  this  introduction  the  differential  equations  of  state  are  now 
presented.     A  nonrotating  frame  with  its  origin  at  the  central  body  is 
assumed. 

s=g(s,  a)  (3-2) 

where 

y 

H-T£  +  a/  (3-3) 

! 

g       (m,a) 
V.   6m  v     s    'J 

and  where  a  is  the  thrust  acceleration,  subsequently  called  the  control 
vector,  \i  is  the  gravitation  constant  for  the  central  body,  and  g     (m,a) 
is  the  mass  rate  equation  for  the  particular  vehicle  under  considera- 
tion.    From  Chapter  II 

2      2 
gm(m,a)  =   -  *-J2-    (VSI)  (3-4) 

2p 

gm(m,a)  =  -   -^L_      (CSI)  (3-5) 

c 

Throughout  the  report  the  variable  m  is  a  normalized  mass,  that  is 
m     =1.     This  obviates  the  necessity  of  selecting  a  value  for  total  ve- 
hicle mass.     The  power  and  the  propellant  flow  are  likewise  normalized 
variables.     That  is 

Exhaust  power  n-fn 

initial  mass 

•        propellant  flow  rate  c\    i\ 

initial  mass 
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3.  5   Linearization  of  the  Equation  of  State 

The  assumptions  of  a  known  reference  trajectory,  known  control 
program  and  of  small  deviations  from  the  reference  trajectory  permit 
the  application  of  linear  perturbation  theory       Thus 


3g  9g 

6s_  =  —      Ss_  +  —      6a 

3s  3a 


From  equation  (3-8)  the  matrices  A  and  B  are  defined. 

°3      h        2 
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The  matrix  G  is  a  symmetric  matrix  of  gravitational  gradients.     The 

31 
extension  to  a  general  gravitational  field  does  not  change  this  property 

However,  numerical  integration  in  such  a  field  requires  an  ephemeris 

for  the  disturbing  bodies.     This  is  a  computational  problem  and  does 

not  affect  the  theory. 

Equations  (3-8)  through  (3  -  15a)treat  the  variations  in  state  velocity 
as  functions  of  the  variations  in  state  and  variations  in  the  control  vec- 
tor, a.     It  is  often  desirable  to  study  the  effects  of  variations  in  per- 
formance parameters  such  as  power,  exhaust  velocity,  propellant  flow, 
etc.     For  this  purpose  one  may  rewrite  the  vector  a  and  the  propellant 
flow  equation  in  terms  of  the  appropriate  parameters.     New  matrices 
A  and  B  will  be  formed  for  this  purpose.     These  cases  are  derived  in 
Appendix  C  since  they  do  not  contribute  to  the  discussion  of  guidance. 
However,  occasionally  in  the  guidance  problem  it  is  desirable  to  work 
with  the  thrust,  f,  as  the  control  instead  of  a.      The  theory  subsequently 
presented  is  applicable  to  this  formulation  also.      For  the  present, 
however,  the  form 

6s  =  A8s  +  B6a  (3-16) 

is  considered  to  be  the  fundamental  formulation  for  the  guidance  prob- 
lem. 

With  the  perturbed  equations  of  state,  as  defined  by  (3-16)  ,  and 
using  the  adjoint  method,  it  is  possible  to  derive  a  state  transition  equa- 
tion. 

3.  6    Method  of  Adjoints 

One  description  for  the  method  of  adjoints  is  obtained  by  consider- 
ing a  set  of  equations  related  to  (3-  16)  by  the  matrix  equation 

A  =  -    A A  (3-17) 

Equation  (3-17)  ,  with  arbitrary  boundary  conditions,  is  said  to  be  ad- 
joint to  (3-  16)  and  the  elements  of  TV  are  the  adjoint  variables.     If  equa- 
tion (3-  16)  is  premultiplied  by  -A.  ,  (3-17)  post-multiplied  by  6j?  and  the 
resulting  equations  are  added,  then  one  obtains 

-^    (  A  6s)     =      AB6a  (3-18) 
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Along  a  given  reference  trajectory,  i.  e     6a  =  O,  equation  (3-18)  has  the 
solution 

A6s_  =  constant  vector  (3-19) 

If  (3-17)  is  integrated  along  the  trajectory,  subject  to  suitable  boundary 
conditions,  then  A  is  a  known  function  of  time  and  the  state  variation 
at  any  time  may  be  determined  from  the  state  variation  at  any  other 
time,  provided  the  reference  control  program  is  used. 

The  adjoint  method  may  be  described  as  the  process  of  introducing 
a  set  of  known  auxiliary  variables,  TV  ,  which  satisfy  (3-17)  and  which 
transform  the  state  variation  at  a  given  time  into  the  state  variation  at 
any  other  time,  provided  a  forcing  function  does  not  exist. 

A  scalar  approach  to  adjoint  equations  and  other  examples  of  their 
application  are  presented  in  Appendix  B. 

A  useful  system  of  equations  closely  related  to  the  adjoint  set  is  a 
set  often  called  the  fundamental  set  or  fundamental  solution.     It  satis- 
fies the  relation 

<£>  =  A  (J)  (3-20) 

Combining  (3-20)  and  (3-17)  by  the  appropriate  pre-and  post-multiplica- 
tion, adding  and  integrating,  yields  the  first  integral 

A  Q  =  constant  matrix  (3-21) 

If  the  boundary  conditions  are  chosen  such  that 

A(tf)    =  I  (3-22) 

$(.0)    -  I  (3-23) 

then  A§>     =     A(0)    =     $(tf)  (3-24) 

This  relationship  is  useful  in  later  discussions 

3.  7    The  State  Transition  Equation 

If  equation  (3-  18)  is  integrated  between  the  times  t1   and  t?J  the 
result  is 

t2 
A2  6s  2  =    A^s^     /      AB6adt  (3-25) 
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Then,  provided  TV  0  is  nonsingular 

*2 


-2  =    A22     A  6-l  +  -^~2l    -f      ^B  5-dt  (3-26) 


h 


-1   . 


A  proof  is  presented  in  Appendix  C  that    _A_       is  not  singular. 

The  product    TV  •       TV     is  called  the  state  transition  matrix  and  is 
denoted  by  T.  .,  where  i  and  j  represent  any  two  times.     Equation  (3-26) 
is  called  the  state  transition  equation.     It  transforms  the  state  varia- 
tion at  t1   and  the  control  perturbation  between  t1   and  t„  into  the  state 
variation  at  t~. 

There  are  three  useful  applications  of  the  state  transition  equation: 

1)  It  may  be  used  as  a  tool  for  studying  the  sensitivity  of  the 
trajectory  to  perturbations  in  launch  conditions  and  to  anoma- 
lies in  engine  performance.     Friedlander       has  performed 
some  investigations  in  this  area.     This  application  is  not 
pursued  further  in  this  report  except  for  derivation  of  appli- 
cable formulations  in  Appendix  C. 

2)  If  engine  performance  is  measured  by  an  accelerometer  or 
other  device  which  can  be  related  to  the  state  equations,  then 
the  state  transition  equation  may  be  used  in  the  navigation  of 
a  spacecraft  by  relating  the  measured  engine  performance  to 
the  state  of  the  vehicle.     This  application  is  discussed  in 
Chapter  IV. 

3)  The  most  important  application  of  equation  (3-26)  is  its  use 
in  guidance.     If,  in  (3-26)  ,  \     is  considered  as  the  final  time, 
then  the  state  and  control  errors  occuring  along  the  trajectory, 
may  be  related  to  the  final  state  variation.     Since  the  mission 
objective  may  be  defined  by  the  final  position  and  velocity  vec- 
tors, (3-26)  meets  the  requirements  for  a  guidance  equation 
provided  that,  for  any  state  variation  6s,  it  is  possible  to  find 
a  solution  a  +  6a  which  will  cause  the  vehicle  to  match  the 
physicalboundary  conditions  at  the  terminal  point.     For  this 
application  the  terminology  "guidance  equation"  is  preferable 
to  "state  transition  equation"  and  will  be  used  to  define  the 
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equation  when  t„  is  the  final  time  and  the  boundary  conditions 
on  the  adjoint  variables  are  chosen  as    A(tf)  =  I.     Thus  the 
guidance  equation  is 


6s     =     A     6s     +    J     AB6adt 
t 


(3-27) 


where  the  initial  time  t  is  any  time  between  the  launch  time 
and  final  time,  at  which  a  new  control  program  is  desired. 

The  integrand  in  (3-26)  and  (3-27)  is  observed  to  be  quite  simple 
when  the  matrices  are  partitioned  and  expanded.     For  the  control  vec- 
tor a 
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Using  (3-10)  and  expanding 
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(3-27a) 


(3-27b) 


where  A  00  is  a  scalar  and  the  remaining    A's  are  three  by  three 
matrices.     The  vectors  are  three  component  null  vectors. 

3.  8   Solution  of  the  FTA  Guidance  Problem 

The  arguments  presented  in  this  section  are  equally  valid  whether 
the  desired  solution  is  the  acceleration  program  a  +  6a  or  a  thrust  pro- 
gram f  +  6f .     A  discussion  of  the  requirements  for  and  the  consequences 
of  interchanging  control  vectors  is  presented  in  Appendix  C.     Since  the 
propellant  expenditure  is  so  easily  written  in  terms  of  the  acceleration 
integral,  especially  for  VSI  vehicles,  it  is  computationally  convenient 
to  work  with  a  as  the  control.     If  a  is  used  in  solving  the  guidance  prob- 
lem, the  set  of  adjoint  functions,    A,  may  be  reduced  immediately  from 
a  seven  by  seven  matrix  to  a  six  by  six  matrix  since  the  equations  of 
motion,  (3-3),  do  not  contain  mass  explicitly  and  the  "mission"  is  to 
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satisfy  terminal  conditions  on  the  equations  of  motion.     However,  when 
working  with  the  CSI  vehicle  it  is  usually  more  convenient  to  use  the 
thrust,  f .     Acceleration  is  used  for  both  vehicles  in  this  report  in  order 
to  be  consistent.     Since  thrust  and  acceleration  are  related  through 
mass,  the  equations  of  motion  using  thrust  as  a  control  must  contain 
the  state  variable,  m.     For  this  reason  the  adjoint  set  cannot  be  re- 
duced a  priori  to  a  six  by  six  matrix.     In  recognition  of  this  fact  and 
since  both  control  vectors  are  useful,  the   A  matrix  is  always  con- 
sidered to  be  a  seven  by  seven  matrix;  with  the  understanding  that  if 
a  is  the  control,  certain  elements  are  identically  zero.     This  does  not 
present  any  problem  in  manipulation  and  need  not  be  considered  fur- 
ther except  when  TV.  and  B  are  to  be  evaluated  numerically.     In  the 
subsequent  discussion,  however,  a  reduced  adjoint  set  is  introduced 
and  is  denoted  by  _A_*.     This  reduction  is  not  associated  with  the  dif- 
ference between  acceleration  and  thrust  formulations.      Both  a  J^   ma- 
trix and  a  _A_  *  matrix  exist  for  each  control  vector. 

From  equation  (3-27)  it  is  clear  that  if  a  state  variation  exists  at 
time  t,  then  in  the  absence  of  control  changes,  a  variation  in  final 
state  will  occur  which  is  given  by 


6sf  =    At  6st  (3-28) 


The  guidance  criterion  is  that  the  final  position  and  velocity  must  satisfy 
mission  requirements.     Presumably  the  reference  trajectory  satisfies 
these  requirements.     Thus  the  position  and  velocity  variations  from  the 
reference  must  be  zero  at  the  final  time  if  FTA  guidance  is  used. 
Stated  mathematically,  it  is  required  that 

f2l 


6sf  =  < 


O 

6mt 


(3-29) 


where  6m„  is  a  small  but  unknown  variation  in  final  mass.     Whatever 
the  value  of  Smf,  it  is  of  no  immediate  concern.     To  reflect  this  it  is 
convenient  to  drop  the  &mf  and  define  the  "miss"    at  the  target  which 
arises  from  the  state  variation  at  time  t  as 

it  ■    A*  6st  (3-30) 
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where  |  is  a  six  component  vector  of  position  and  velocity  variations 
and  TV*  is  a  six  by  seven  matrix  obtained  by  deleting  the  seventh  row 
of  A  .     Clearly  6mf  can  be  obtained,  once  the  corrected  control  is 
known,  from 

h 

6m    =AJ+6S++     f     A  J   B  6a  dt  (3-31) 

I  (X      —  t         ,J       —    (  — 

where  A  7  denotes  the  deleted  seventh  row  of  the  adjoint  set. 

Since  the  state  error  6s,  is  assumed  to  be  small,  so  that  linear 
theory  is  valid,  it  follows  that  a  small  correction  to  the  reference  con- 
trol program  will  be  sufficient  to  null  the  miss  vector  in  the  remaining 
flight  time.     The  foregoing  is  true  provided  the  vehicle  has  sufficient 
thrust.     That  is,  it  is  not  "in  extremis".     Assume  it  is  not.      Then  the 
miss  may  be  reduced  at  a  rate  such  that 

|  =    A  *  B  6a  (3-32) 


or  by  integrating,  such  that 

O  -  lt  =     J     A*  B  6a  dt  (3-33) 


*f 


That  the  control  6a  is  not  unique,  except  for  a  |_    which  requires 
maximum  thrust  continuously,  is  evidenced  by  the  fact  that  for  ballistic 
guidance,  ideally,  only  two  corrections  are  needed  to  null  position  and 
velocity  error;  a  midcourse  correction  to  correct  position  and  a  termi- 
nal correction  for  velocity.  Thus  for  continuous  thrusting  vehicles  an 
infinity  of  solutions  exists.     The  criterion  for  selecting  a  unique  con- 
trol has  already  been  given;  the  acceleration  integral  J  or  J  *  must  be 
minimum.     Mathematically,  the  problem  for  VSI  vehicles  may  be 
stated:     It  is  required  that 

Tnl  tf 

J—  I      =  £      +      /       A¥B  6a  dt  (3-34) 

\°)  t 

and  that  J  is  a  minimum,  where 

J  =  J'    <a+6a)T(a+6a)     flt  (3_35) 

t  Z 
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For  convenience  define  the  optimal  control,  a 

a°  =  a  +  6a  (3-36) 

For  a  VSI  vehicle  with  no  constraints  on  a     the  set  (3-34)  and  (3-35) 
can  be  solved  by  a  direct  application  of  the  calculus  of  variations.     The 
vector  space  around  the  reference  trajectory  is  explicitly  assumed  to 
be  flat  to  first  order.     Consequently  TV  and  B  are  invariant  between 
neighboring  trajectories.     Thus 

f  tf  1\ 

dt  +  I       <\   1  +    I      -A-*  B(a°  -  a)  dt   ;  )     =0 

it  J/  (3-37) 

where  n  is  a  vector  of  constant  Lagrange  multipliers  and  the  subscript 

on  |_    is  dropped.     Expanding  (3-37)  one  obtains 

T 
a°      =  -  £     A*  B  (3-38) 

The  vector  ir  is  eliminated  with  help  of  equation  (3-33). 
r     .      ~~  tf 

<-J.      =  J_  -     /  A*  b{bT  A*T  £  +  a  I     dt  (3-39) 

|oJ  t  I  / 

-  7T  =  M"1  (77  -  £  )  (3-40) 


where 


tf 


Thus 


M=     /    A-BB1    A*      dt  (3-41) 

t 

T]  =     /   A*  B  a  dt  (3-42) 

t 


:  b1    A-1   M"1    (77  -Jj  (3-43) 


The  solution  (3-43)  is  valid  provided  the  reference  trajectory  is  suffi- 
ciently close  to  an  optimal  that  6a  is  small.     It  is  unique  provided  M  is 
not  singular.     A  proof  is  presented  in  Appendix  C  for  the  existence  of 
M"1. 

A  physical  interpretation  of  the  quantities  M  and  77,  defined  by 
(3-41)  and  (3-42),  is  valuable  in  understanding  the  solution.      First  ob- 
serve that  I  has  two  interpretations.     From  (3-30),  £  is  the  miss 
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arising  from  a  state  variation,  6s  ,,  when  the  reference  control  is  used. 
From  (3-33),   -  £  may  be  interpreted  as  the  miss  that  results  from  a 
control  variation,  6a,  occuring  after  time  t  but  with  6s,   =  £.     Similarly 
from  (3-42),   -  77   is  interpreted  as  the  particular  miss  that  arises  if 
6a  =  -a;  that  is,  coast  is  initiated.     If  the  miss,  £_,  due  to  a  state  varia- 
tion and  the  miss,   -  77,  due  to  initiating  coast  are  such  that  (l_-7?  )    =  £, 
the  optimal  control  is  the  null  vector.     Although  such  a  solution  will 
generally  violate  linearity  assumptions,  it  will  indeed  result  in  a  mini- 
mum J,  namely  zero. 

If  the  product  _A_  *  B  is  interpreted  as  the  sensitivity  of  the  final 
miss  to  a  unit  impulse  in  each  component  of  the  control  at  time  t,  then 
M  may  be  considered  as  a  weighted  total  sensitivity  of  the  miss  to  a 
unit  control  applied  continuously  between  t  and  tf.     In  a  physical  sense 
the  procedure  computes  the  vector  sum  of  the  miss  using  the  reference 
control  vector  and  the  miss  using  a  null  control  vector;  then  selects 
the  control  at  each  point  according  to  the  sensitivity  of  the  miss  at  that 
point. 

It  is  interesting  to  note  the  solution  which  results  if,  instead  of  J, 

the  minimization  criterion  is 

tf    c    T, 
1     6a      6a 

S  =     J       — =  dt  (3-44) 

t  2 

Proceeding  as  in  equations  (3-37)  through  (3-41)  but  solving  for  6a,  the 
result  is 

a°  =  a  +  6a  =  a  -  BT   A  *T  M'1  J_  (3-45) 

T  T       - 1 

Comparing  (3-45)  with  (3-43)  it  is  observed  that  for  a  =  B     W*      M     77 

the  results  are  the  same.     The  preceding  expression  for  a  satisfies 

(3-42)  therefore  for  small  variations  around  the  reference  trajectory 

the  two  methods  result  in  the  same  control.      This  is  only  true  for  VSI 

trajectories,  however. 

Before  proceeding  to  consider  CSI  solutions  it  is  desirable  to 
examine  the  consequences  of  the  reference  trajectory  being  optimal 
for  the  initial  launch  point      If  the  reference  is  initially  optimal  but  a 
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state  variation  occurs  such  that  a  miss,  £,  will  result  from  use  of  the 
reference  control  program,  the  reference  control  ceases  to  be  an 
optimum  because  the  boundary  conditions  are  not  satisfied.     Also  the 
particular  formulation  used  for  the  problem  does  not  result  in  a  singu- 
lar M  matrix  as  is  often  the  case.     Therefore,  in  this  study,  the  only 
consequence  of  using  an  optimal  reference  is  the  assurance  that  6a 
will  not  violate  linearity  assumptions. 

A  new  approach  is  now  presented  which  can  be  used  for  both  VSI 
and  CSI  control.  To  illustrate  this  method  the  VSI  problem  is  solved 
again  but  with  a  constraint. 

Assume  that  thrust  is  not  to  exceed  an  amount  f   .      Therefore 

o 

thrust  acceleration  cannot  exceed  f  /mm    .     To  be  consistent  with  the 

o'  o 

normalized  mass  variable,  m,  used  in  this  report  denote  f  /m     as  the 

^  o'      o 

initial  acceleration  limit  a   .     The  constraint  may  now  be  written  as 

o  J 

a 

a<  —  (3-46) 

—  m  v  / 

1  R 

Following  Kalman       and  using  Pontryagin's  principle  (Appendix  D), 

form  the  scalar  Hamiltonian,  H 

T 
o        o 
a        a  q-, 

H  =  ~ —    +    ul  |  (3-47) 

2 

The  theory  states  that  if  the  cost,/((a    )   / 2)dtis  to  be  a  minimum,  then  for 
each  point  along  the  trajectory  H  must  be  a  minimum,  where  v_  is  an 
unknown  vector,  often  called  the  costate,  which  satisfies  the  adjoint 
relationship,  and  _£    is  the  state  velocity.     The  linearization  in  this  chap- 
ter permits  considerable  simplification.     Using  (3-32),  £_  is  interpreted 
as  a  variable  representing  the  velocity  of  the  final  state.     The  variable 

« 

|_,  is  the  initial  condition  for  £    and  is  a  function  of  the  lower  limit  of 
integration,  t.     J,  is  evaluated  from  the  vehicle  state,   6s,,  using  (3-30). 
Because  £_  is  a  variable  in  state  space  at  the  final  point,  v_  is  a  constant 
vector.     In  particular  it  is  the  final  value  of  the  general  time  varying 
costate  vector. 


Therefore,  using  (3-36) 


\    =    A*  B  (a°  -  a)  (3-48) 
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Further  T 

o        o 

cL  3. 

H  =— —    +  v1   A*  B  (a°  -  a)  (3-49) 

2 

T 
o        o 

cL  3. 

H  -  -  vT    A*  Ba  +— —    +  v^    A*  Ba°  (3-50) 

2 

The  first  term  in  (3-50)  is  a  constant  independent  of  a   ,  thus  H  is  a 
minimum  when  the  sum  of  the  last  two  terms,  H..,  is  a  minimum. 

The  control  a     is  to  be  determined  such  that  H.  is  a  minimum. 

T 
o        o 
a        a  ,-p 

minimize  H1   =  +  y     A*  Ba  (3-51) 

2 

a 

subject  to  a  <     (3-46) 

—      m 

now  let  ^T  A*  B  =  qT  (3-52) 

Since  q  and  a     are  both  three  component  vectors,  one  may  be  obtained 
from  the  other  by  a  scalar  multiplication  and  a  rotation.     Thus 

a°  =  y  C  q  (3-53) 

where  y  is  a  positive  scalar  and  C  a  coordinate  rotation.     Inserting 
(3-53)    into  (3-51) 

2     2 
H     =2— ^ —  +yq      Cq  (3-54) 

2 

T 
For  any  value  of  y  ,  H     is  minimum  if  q     Cq  has  maximum  magnitude 


and  is  negative.     But 


qTCq|    <     |q|    |Cq|  (3-55) 


where  equality  holds  for  C  =  I.     It  is  apparent  that  H     is  a  minimum 

2 
only  if  C  =  -  I  and  (y    /2  -  y)   is  a  minimum. 

Hx  =  P—   -  y)     q2  (3-56) 
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a 
o            o 

q 

q> 

a 
o 

m 

q 

m 

o 
a     =  -  q 

q£ 

a 
o 

m 

y  -  min 

f. 

a 
o 

) 

Thus  when  constraint  (3-46)  is  applied,  H     is  a  minimum  if  and  only  if 

(3-57) 

(3-58) 

(3-59) 
mq 

or 

a°  =   -  y  BT   A*T  v_  (3-60) 

Comparing  this  solution  with  the  variational  solution,  equation  (3-38), 
one  observes  that  except  for  the  constraint  y  the  solutions  are  identical 
and  v_  -  it.     y  may  be  interpreted  as  a  switching  function  which  carries 
the  thrust  restriction  (3-46). 

Having  gained  confidence  in  the  Hamiltonian,  now  apply  the  method 
to  the  CSI  transfer  for  which  the  classical  approach  is  degenerate.     The 
appropriate  cost  is  t^  o 

J*  =  /    — dt,  (2-19) 

t         2m 

Minimize  o 

a     a 

H  =  — +  ^T    A*  B  (a°  -  a)  (3-61) 

2m 

a 

subject  to    a  <   (3-46) 

m 

T 
Again  v    J\J?  B  a  is  a  constant,  therefore 

■    •      •        tt  o       o   .      T     o  ,0    cn. 

minimize  H     = a     +  q      a  (3-62) 

2m 
Where  q  is  defined  as  before. 

The  previous  argument  holds,  to  yield 

a     =  -  y  q  (3-63) 


and 

H 


r-r(^--q)q  (3-64, 


40 


In  this  case  H     is  a  minimum  if  and  only  if 


where 


o 

a     = 


y  =  0 


y  BTA*T  v 


(3-65) 


y 


o 


mq 


a 


q 


< 


> 


o 


2m 
a 


o 


2m 


(3-66) 


For  the  case  q  =  a   /2m,    y  is  actually  indeterminate.     However  this  is 
not  an  important  consideration  for  this  thesis  since  equality  holds  only 
for  infinitesimal  time  periods. 

Equations  (3-65)  and  (3-66)  yield  the  "bang-bang"  solution  char- 
acteristic of  optimal  trajectories  for  which  the  cost  is  a  linear  function 
of  the  control.     Again,  the  solution  is  strictly  valid  only  for  perturba- 
tions around  a  reference  trajectory.     In  order  to  complete  the  CSI  solu- 
tion it  is  necessary  to  evaluate  the  constant  vector  v_.     The  evaluation 
is  more  difficult  in  the  CSI  case  than  for  VSI  because  the  control  pro- 
gram is  discontinuous,      v  is  evaluated  by  first  applying  the  boundary 
conditions  (3-34)  and  using  (3-36). 

t. 


'f 
0  =  i+     J    A*  B  (a°  -  a)  dt 

t 


(3-67) 


Assume  that  a     and  a  differ  only  because  their  respective  values  of  y 
and  V_  are  different.     The  assumption  merely  implies  that    A*  B  is  the 
same  for  neighboring  trajectories.     With  this  assumption  rewrite 
(3-65)  as 

O  O    DT  .,,,T    „0  /q     e0\ 

a     =  -  y     B   A*      v_  (3-68) 

for  the  corrected  control  program,  and  let 

a  =  -  y  BTAi;T  v  (3-69) 

represent  the  reference  control  program.     Then  (3-67)  becomes 

O  =  |  - '  /  A*  b  btA*t  (r°  _^°  -  yv  )  dt  (3-70) 

t 
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which  may  be  rewritten  as 


|=     /     A*BBT    A*        A(y^)dt 
t 


(3-71) 


The  increment  A  (y  ^)  represents  the  difference  in  the  control  pro- 
gram between  the  reference  trajectory  and  the  corrected  trajectory. 

Since  y  is  a  function  of^,  a  solution  may  be  obtained  by  formally 
differentiating  and  solving  (3-72)  for  6.  \j. 


d 


9  v     t 


/       A*  B  BT  A*Ty  u  dt 


t±v 


(3-72) 


The  function y   is  discontinuous  at  switch  points,  which  occur  at  times 
t,    in  the  interval  tf  -  t.     Consequently  the  integral  (3-72)  must  be 
separated  into  regions  of  coasting  and  thrusting  and  Leibniz's  rule 
used  for  the  differentiation.     A  term  of  the  form 


±  A*  b  bt  A*Ty  v 


at, 


Lk  a, 


A  v     results  for  each  switch  point,  where 


plus  is  for  switch  off  and  minus  is  for  switch  on. 

The  term  dy  /d  y_  is  evaluated  in  the  continuous  regions  from  the 
definitions  (3-52)  and  (3-66).  The  term  8t, /3_^is  evaluated  by  con- 
sidering 


q  =  q  (t,  V  ) 


From  (3-66 


q  (tk,^)   = 


o 


2m 


(3-73) 


(3-74) 


Differentiating  (3-74)  with  respect  toj^,  one  obtains 

■    M_    =  o 


3q  k 


at         dv_ 

a  v 


a^ 

3q 
du 

at 


(3-75) 


(3-76) 
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The  derivatives  of  y   and  q  are  evaluated  in  Chapter  VI  when  it  becomes 
necessary  to  specify  the  explicit  form  for  computation. 

Solving  for  \v ,  one  obtains 


A  v  =  M*"1  £  (3-77) 


where 


yl  +  v  — - 


M*  =  /A*  B  BTA*T 

t 


a  yj 


dt  ±A*  bb  A     yv 


The  constant  vector  y     is  obtained  from 

v     -  v  +  &P 


where 


atk 

*k 

dv 

(3- 

-78) 

(3- 

-79) 

and 


v  =  -  M"%  (3-80) 

M  =  /   y  A*  B  BT  A*T  dt  (3-81) 

t 

Z7  =  /  A*  Ba  dt  (3-82) 

t 

The  values  for  y     are  obtained  from  equation  (3-66)  by  using  the  com- 
puted v_    in  q. 

With  all  quantities  defined  in  the  preceding  equations,  the  corrected 
control  program  is  given  by 

a°  =  y°  BTA*T  (m'%  -  M*_1i)  (3-83) 

3.  9   Application  of  the  Guidance  Theory 

In  order  to  use  the  theory  of  section  3.  8  for  the  guidance  of  space 
vehicles  it  is  necessary  to  compute  the  quantities  in  equation  (3-43) 
for  VSI  guidance  or  in  equation  (3-83)  for  CSI  guidance.     Since  (3-43) 
may  be  regarded  as  a  special  case  of  (3-83),  only  the  latter  is  discussed. 

If  the  reference  trajectory  is  known,  then  the  elements  of  the  ma- 
trices B,  A,  M,  and  M*  and  the  components  of  the  vector  r\   can  be 
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computed  for  any  point  along  the  trajectory      Whether  these  quantities 
are  precomputed  and  stored  prior  to  launch  or  are  computed  on  board 
as  needed  will  be  determined  by  the  state  of  computer  technology  when 
the  vehicle  is  designed.     This  problem  is  not  relevant  to  the  present 
discussion.     In  either  case,  the  first  step  in  computing  the  corrected 
control  is  to  establish  that  a  state  variation,  6s  ,,  exists.     Methods  for 
processing  measurement  data  for  this  purpose  are  discussed  in  Chap- 
ter IV.     Then  the  final  state  error  is  determined  from  (3-30). 


I  -   A*  6s  t 


(3-30) 


The  vector  y_     is  computed  from  the  known  quantities  and  the  final  state 
error  using  (3-77)  through  (3-82). 


iP  =  M*"1^   -  M"1  T2 

o  o 

From  (3-52)  and  (3-66)  q     and  y     may  be  computed. 


o 


BTA- 


<^T  u° 


y 


o 


y 


o 


o 


k 


mq 


o 


q°< 


q°> 


o 


2m 
a 


(3-84) 

(3-52) 
(3-66) 


o 


2m 


Then  from  (3-68) 


o 

a     = 


}' 


°BiA- 


*T     o 


(3-68) 


This  result  may  be  programmed  into  the  vehicle  control  system  to 
implement  FTA  guidance. 

3. 10  Discussion  of  the  VTA  Guidance  Problem 

In  Chapter  II  it  is  shown  that  in  field-free  space,  the  minimum 

2       3 
acceleration  integral  for  VSI  vehicles  is  J  -  6L  /T   .     A  VTA  guidance 

scheme  which  minimizes  J  for  a  given  L  in  field-free  space  will  there- 
fore select  an  infinite  transfer  time  unless  constrained.     If  the  planets 
were  all  in  coplanar,  circular  orbits  and  the  reference  trajectory  lay 
in  this  plane,  a  similar  result  would  be  obtained  in  the  solar  system. 
Because  of  the  inclination  and  ellipticity  of  planetary  orbits,  local  min- 
ima   of     J     will  occur  which  depend  upon    planet    and    spacecraft 
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orientation.     A  complete  relaxation  of  the  time  constraint,  tfJ  in  the 
theory  of  sections  3,  8  and  3.  9  would  require  that  the  new  arrival  time 
correspond  at  least  to  a  local  minimum  of  J.     Therefore,  unless  tf  is 
approximately  equal  to  this  propellant-optimal  arrival  time,  linear 
theory  cannot  be  used  to  obtain  propellant-optimal  VTA  guidance. 

However,  it  is  possible  to  restrict  the  change  in  arrival  time,  use 
linear  theory  and  obtain  useful  solutions.     A  case  where  this  is  impor- 
tant occurs  for  a  vehicle  "in  extremis".      That  is,  insufficient  power 
(or  thrust)  is  available  to  null  |  at  time  t.~.     The  result  obtained  is 
dependent  upon  the  method  of  restricting  At  and  upon  any  simplifying 
assumptions.     Consider  the  following  example  which  is  based  upon  the 
assumptions: 

1)  The  target  point  relative  to  the  planet  center  is  unchanged. 

2)  The  final  velocity  relative  to  the  target  point  is  unchanged. 

3)  The  thrust  acceleration  af  =  a  (tf)  is  constant  over  the  interval. 
The  motion  and  position  of  the  target  point  at  the  time  tf  +  At  are 

r  T  (tf  +  At)  =  r  T  (tf)    +  vp  At  (3-85) 

vT  (tf  +  At)  =  vT  (tf)  =  yp  +  vR  +  gT  At  (3-86) 

where  v  and  vR  are  the  planetary  motion  and  the  desired  relative 
velocity,  respectively,  and  g  T  is  the  solar  gravity  at  the  planetary 
radius. 

For  the  mission  to  be  accomplished,  the  vehicle  position  and  veloc- 
ity   must  equal  r  „  and  v  T  at  time  t~  +  At       The  vehicle  position  and 
velocity  in  terms  of  the  reference  trajectory  are 

r  (tf  +  At)  =  r    (tf)  +  6r  (tf)  +  v  (tf)    At  +  af     (~  \     +  Sv  (tf)  At 

(3-87) 
y  (tf  +  At)  =  y  (tf)    +  6v  (tf)    +af(At)    +  gT  At  (3-88) 

Solving  equations  (3-85)  through  (3-88)  for  6rf  and  6y_  one  obtains 
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6£f  =  _  XR   At  ~  6Xf  ^~  (3-89) 

6vf  =  -  af  At  (3-90) 

For  the  linear  approximation  assume  6v„/2   «  vR.     Equations  (3-89) 
and  (3-90)    may  be  inserted  in  the  guidance  equation. 

-<(      ">At     =  £  +    J   A*  B  6a  dt  (3-91) 

UfJ         * 

Define  a  new  miss  vector  |>;c  such  that 


£*  =   <  >     =<  >      +<  >    ^  (3-92) 


Then  by  selecting  At  such  that  the  uncorrected  final  position  error, 

!*       ,  is  a  minimum,  the  vehicle  will  attain  the  terminus  but  at  the 
time  tf  +  At  and  with  the  least  departure  from  the  actual  trajectory. 
The  result  is 

T  t 
At  =  ~ — -  (3-93) 

VR 

Equation  (3-93)  may  be  inserted  into  (3-92)  and  the  optimal  control 
program  found  using  the  FTA  procedures  of  sections  3.  8  and  3.9. 

If  the  target  point  is  the  planetary  sphere  of  influence  and  the  de- 
sired relative  velocity  is  zero,  obviously  the  above  solution  is  invalid. 
Other  assumptions  may  be  used  to  treat  such  cases.     The  assumptions 
to  be  used  and  the  criterion  for  restricting  At  may  be  changed  to  suit 
the  purposes  of  the  investigator.     In  general,  the  propellant -optimal 
VTA  problem  is  not  readily  treated  by  linear  methods  unless  additional 
constraints  are  used.     Exploring  the  numerous  possibilities  that  arise 
will  be  left  for  future  investigations. 


■Hi 


CHAPTER  IV 


ESTIMATE  OF  THE  STATE  VECTOR 


4.  1   Introduction 

In  this  chapter  a  method  of  determining  a  "best  estimate"  of  the 
vehicle  state  vector  is  presented.  The  term  "best  estimate"  is  used 
to  designate  the  estimate  for  the  state  vector  computed  by  processing 

redundant  measurement  data  through  a  statistically  optimum  filter, 

32 
Potter  and  Stern   '    have  shown  that  all  of  the  three  commonly  used 

methods  of  processing  redundant  data,,  namely    maximum  likelihood, 
minimum  error  ellipsoid,  and  minimizing  a  characteristic  scalar 
parameter  all  result  in  the  same  filter.     They  have  shown  further  that 
for  each  unbiased  optimum  filter,  there  exists  an  associated  biased 
optimum  filter  which  produces  a  smaller  error  ellipsoid  for  the  esti- 
mate 

Since  all  of  the  methods  result  in  the  same  filter,,  in  this  thesis 
the  method  which  presents  the  least  mathematical  complexity  is  used. 
The  method  involves  a  slight  variation  of  the  minimum  error  ellipsoid 
technique. 

4.  2    General  Remarks 


To  establish  a  corrective  thrust  program  it  is  desirable  to  have 
accurate  knowledge  of  the  state  vector  at  the  time  the  program  is  to  be 
initiated.     It  is  clear  that  given  perfect  knowledge  of  some  prior  state 
vector  and  of  vehicle  performance  the  prediction  problem  reduces  sim- 
ply to  application  of  the  state  transition  equation  between  the  time  for 
which  the  state  is  known  and  the  time  for  which  the  prediction  is  de- 
sired.    It  is  equally  clear  that  perfect  knowledge  of  a  process  seldom 
exists.     Thus  the  problem  becomes  that  of  using  existing  information 


An  unbiased  filter  will  produce  a  true  value  for  the  state  vector  if 
all  measurements  are  free  of  error.     A  biased  filter  is  biased  in  favor 
of  the  a  priori  or  prelaunch  expectation  of  the  second  moment  of  the 
state  vector  probability  density, 
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in  such  a  way  that  its  inaccuracies  have  minimum  effect  on  the  guidance 
decisions. 

For  ballistic  vehicles  the  determination  of  position  at  two  points 

along  the  trajectory  suffices  to  specify  the  entire  trajectory      Several 
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investigators  including  Battin       and  Stern      ,  have  studied  the  effects 

of  using  additional  fixes  and  of  selecting  the  geometry  of  individual 
celestial  measurements  to  reduce  the  error  in  the  trajectory  determina- 
tion.    The  techniques  developed  in  those  studies  may  be  extended,  in 
some  cases  without  change,  to  the  present  work.     The  primary  differ- 
ence between  the  ballistic  and  continuous -thrust  vehicles,  is  that  the 
future  portion  of  the  actual  continuous-thrust  trajectory  can  never  be 
completely  determined  on  the  basis  of  its  history  alone.     This  fact  is 
due  to  the  possible  occurrence  of  random  changes  in  the  thrust  vector. 
Another  difference,  which  is  more  easily  treated,  is  the  first  order 
dependence  of  the  trajectory  upon  vehicle  mass.     This  dependence  may 
be  handled  by  measuring  propellant  state  at  discrete  intervals  as  well 
as  making  celestial  measurements 

Thus  for  the  continuous-thrust  spacecraft,  state  determination  from 
on-board  measurements  may  be  separated  into  two  related  problems: 

1.  determining  the  state  history  by  suitable  filtering  of  all  meas- 
urements that  have  been  made,  and 

2.  predicting  future  values  of  the  state  vector. 

The  second  of  these  is  dependent  upon  but  is  not  uniquely  determined 
by  the  first. 

The  use  of  celestial  sightings,  as  in  the  case  of  ballistic  transfer, 
is  sufficient  to  determine  the  spacecraft  trajectory.     However  such 
measurements  are  not  sufficient  to  specify  the  entire  state  history  nor 
the  variational  history  of  the  control  vector  which  contributed  to  the 
state  change.     Because  of  this,  prediction  of  the  state  vector  solely 
from  periodic  celestial  sightings  and  propellant  measurement  will  be 
subject  to  larger  uncertainties  than  if  the  prediction  includes  meas- 
urements of  the  control  vector  as  well. 
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4.  3    Unbiased  Estimate  of  Present  State 

The  term  "present". .state  will  be  used  to  denote  the  state  vector 
computed  without  regard  to  a  computational  time  lag,  thus,,  it  is  the 
state  at  the  time  of  the  most  recent  measurements. 

The  measurements  which  will  be  processed  for  the  estimate  are 
1)  periodic  determination  of  position,  2)  propellant  state,  and  3)  con- 
tinuous measurement  of  engine  performance.     It  appears  quite  feasible 
to  derive  a  method  of  including  raw  celestial  observations  in  the  compu- 
tation process  without  explicitly  deriving  the  position  vector.     How- 
ever, this  would  contribute  nothing  to  the  present  argument  and  will 
be  left  as  a  subject  for  future  study      Consequently,  the  computed  com- 
ponents of  position  variation  will  be  considered  as  a  "measurement". 
Since  velocity  is  impractical  to  measure  directly,  except  near  planets, 
the  velocity  will  be  derived  from  successive  position  measurements. 

The  time  of  the  present  state  estimate  will  be  chosen  as  occuring 
between  celestial  "fixes".     If  the  state  at  the  time  of  a  fix  is  desired, 
such  an  estimate  will  be  a  special  case  of  the  more  general  problem. 
This  approach  is  justified  on  the  basis  that  celestial  observations  may 
be  separated  by  several  days  but  current  engine  data  is  always  avail- 
able. 

In  the  following  development  measured  quantities  will  be  denoted 
by  the  "tilde"  (~)  and  the  estimated  quantity  by  a  carat  (A)- 

The  relation  between  the  measurement  quantities  and  the  state  may 
be  written  as: 

N6s  =  6q  (4-1) 

where       N       is  a  k  by  seven  deterministic  matrix  relating  6s  to  6q, 
6  s     is  the  state  to  be  computed 

6q     is  a  column  vector  whose  k  components  are  the  meas- 
urement data 

The  matrix  N  may  be  derived  directly  from  the  state  transition 
equation  written  between  the  time  for  which  the  estimate  is  desired, 
t  =  t     and  any  previous  time    t  =  t. 
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TV    •     For  simplicity,  the  integral  term  in  equation 
(4-2)  may  be  considered  as  a  vector. 
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and  the  transition  matrix  T.     as  a  partitioned  matrix  of  partial  deriva- 
tives.    Thus,  if  only  the  position  vector  at  t.  is  considered,  one  obtains 
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where  the  matrix  of  partial  derivatives  consists  of  the  first  three  rows 
of  T.    .     Equation  (4-4)  may  be  written  for  any  number  of  times 

(i  =  n  -  1,  n  -  2  .  .  .'.  ),     An  analogous  equation  may  be  written  for 
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the  propellant  measurement  and  the  results  arranged  in  the  form  of 
equation  (4-1).     One  might  expect  the  propellant  measurement  to  be 
less  critical  and  less  subject  to  error  than  the  other  measurements. 
Consequently  it  may  be  necessary  to  include  only  a  few  measurements 
of  propellant,  perhaps  one  or  two. 

Let  us  consider  now  the  problem  of  finding  an  unbiased  filter,  F„, 
which  minimizes  the  error  in  the  state  vector  and  satisfies  the  rela- 


tion 


6g=   FQ6q 


(4-5) 


If  the  uncertainty  in  the  state  vector  is  u  and  the  error  vector  associat- 
ed   with  measurement  is  e ,  then 


=  6 


N  6s  +  Nu 


(4-6) 


Since  the  true  measurements  and  the  actual  state  vector  satisfy  equa- 
tion (4-1).  then 


e  =  Nu 


(4-7) 
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The  equation  for  the  error  ellipsoid  associated  with  the  measure- 
ment vector  is  given  by 


where 


eT  E"1  e    =  1  (4-8) 


E  =  <  e    €T>  (4-9) 


is  a  k  by  k  matrix  and  the  brackets  denote  the  mean  value..     Equations 
(4-6),  (4-7),  and  (4-8)  may  be  combined  to  yield 

u1   N1   E"1  Nu  =  (6q  x    -  6  s1   N1)  E_i  (6q  -  N6  s)        (4.10) 

If  the  partial  derivative  of  equation  (4-10)  with  respect  to  the  compon- 
ents of  the  state  vector  are  set  equal  to  zero,  one  obtains: 

■ 

fd-    \        T      -1  T      -1       ~ 

2      — -        N1   E        Nu  -   -   2  N1   E        (6q  -  N6s)  -  0 

The  solution  to  equation  (4-11)  is  the  estimate,  6  s. 


(4-11) 


6s  =  (NTE-1N)~     NT  E"1  6q  (4-12) 

Thus  for  the  unbiased  filter,,  equation  (4-12)  yields: 

T      -1        ~1      T      -1 
FQ  =  (N      E        N)        N      E  (4-13) 

It  is  not  surprising  that  the  filter  of  equation  (4-13)  is  the  same  filter 
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derived  by  Potter  and  Stern  using  the  method  of  maximum  likelihood 

4.  4   Biased  Estimate  of  Present  State 


The  proof  that  a  biased  estimate  will  result  in  a  smaller  ellipsoid 

of  error  will  be  omitted  since  this  topic  is  well  covered  in  the  litera- 
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ture  In  this  section  only  the  method  of  obtaining  the  biased 

filter  from  the  unbiased  filter  will  be  presented. 

The  optimum  biased  filter,  F^RJ  is  obtained  by  deleting  the  last 
seven  columns  from  an  associated  unbiased  filter  which  is  computed 
using  a  fictitious    measurement  of  present  state.     Thus  the  measure- 
ment error  vector  to  use  for  computation  of  the  error  covariance 
matrix  E  is 
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v    -  Ss  ^ 

The  last  seven  rows  of  the  N  matrix  must  then  be  the  seven  by  seven 
identity  matrix.     The  resulting  unbiased  filter  F„  computed  from 
equation  (4-13)  but  with  the  last  seven  columns  deleted  defines  the 
optimum  biased  filter,  FnR 

4.  5   Covariance  Matrix  of  Measurement  Error 


The  matrix  E  is  extremely  important  in  the  computation  of  the 
state  estimate.     Since  it  involves  the  measurement  errors  from  both 
discrete    and  continuous  measurements,  it  is  worthy  of  closer  examina- 
tion. 

In  section  4.  3  a    subvector    of  the  measurement  vector  Sq  was 
written  in  the  form 


6ai 


6  r  .  +  A  r  . 
—  l  — l 


From  equation  (4-3)  it  is  apparent  that 
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£.  =  [  I3  Os  03]    Aj         /    AB  6a  dt 


(4-15a) 


(4-15) 


where  I„  is  the  three  by  three  identity,  0„  is  a  three  by  three  null 
matrix  and  O  „  a  three  component  null  vector.      Thus  a  measurement 
of  position  consists  of  two  parts:     the  celestial  fix,   &£ .  ,  at  time  t  -  t. , 
and  an  additional  vector  Ar  .  which  contains  the  integrated  engine  vari- 
ations from  t.  to  t    .     As  a  consequence,  the  error  vector  at  time  t. 
in  ^  i 

will  also  consist  of  two  parts, 

'        ti 

(4-16) 
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which  for  simplicity  will  be  written  as 


52 


*n 


.  =  e    .  +    f    D.  €      dt  (4-17) 

1     —  ri     ,J        l—  a 


4i 

where  e     .  is  the  error  in  determining  position  and  €      is  the  error  in 
—  ri  te  ^  —a 

measuring  the  engine  quantity. 

It  appears  justifiable  to  assume  that  measurement  errors  of  dif- 
ferent types,  i.  e.    position,  propellant  flow,  acceleration,  are  indepen- 
dent and  have  zero  means.     Further,  that  the  position  measurement 
errors  for  two  different  times  are  independent.     The  assumption  of  a 
zero  mean  for  measurement  errors  does  not  result  in  a  loss  of  gener- 
ality, however  it  greatly  simplifies  the  mathematics 

With  the  above  assumptions,  the  E  matrix  includes  terms  of  the 
following  form: 

1)      On  the  diagonal 


(4-18) 
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2)      Off  the  diagonal 

r11  rn     T      T 

<     /     D.  e     dt     /      e       D.     dt  >  (4-19) 

,J        l  —  a        ,J     —  a       i 

o  o 

Using  the  theory  for  handling  random  processes         the  integral 


terms  are  easily  reduced  to -the  form 

t  t 

n  n 

J     D.  (t    )    J     6   (t  t2)  D1   (t2)  dtx  dt2  (4-20) 

t  •  t .  J 


where  6   (t-tp)  is  a  diagonal  matrix  of  autocorrelation   functions.     To 
proceed  it  is  necessary  to  make  some  assumption  about  6  .     The  most 
easily  justified  assumption  is  that  the  engine  measurement  errors 
each  contains  much  higher  frequency  components,  than  does  the  matrix 
D.  and  each  is  uncorrelated  except  over  short  intervals       Thus,  over 
the  period  of  integration  the  measurement  errors  approach  a  "white 
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noise"  distribution.     With  this  assumption  9  is  reduced  to  a  diagonal 
matrix  of  constants  representing  the  measurement  error  variance 
multiplied  by  the  Dirac  delta  function. 

r2 
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6(VV 


(4-21) 


(Measurement  of  propellant  flow  and  exhaust  velocity  are  used  as  ex- 
amples in  (4-21).)     Designate,  for  simplicity 

o  =  E  6(t2-ti)  {4_22) 

where  2j  is  a  diagonal  matrix  of  measurement  variance. 

One  integration  of  the  terms  in  (4-18)  and  (4-19)  may  be  performed 
to  yield 

f    D.   Y,   DT  dt  (4-23) 

where  the  interval  k  to  n  represents  the  shorter  of  the  intervals  i  to  n 
or  j  to  n.      The  E  matrix  formed  from  these  elements  is  a  symmetric 
positive  definite  matrix. 

The  necessity  for  including  engine  measurements  in  the  estimate 
of  state  for  low-thrust  vehicles  is  made  evident  by  equation  (4-24) 
Note  that  the  measurement  error  vector,  when  position  data  alone  is 
used,  would  contain  integrals  of  the  actual  engine  variations.     That  is 

m 

6.  =  €    .  +     f    D.    6a  dt  (4-24) 

—  l     —  ri        J         i     — 


Assuming  that  the  engine  measurement  errors  are  much  smaller  than 
the  engine  variations  to  be  measured,  equation  (4-24)  would  result  in  a 
sizeable  increase  in  the  error  ellipsoid.     In  addition,  as  the  interval 
i  to  n  increases,  the  value  of  the  celestial  fix  at  time  t.  is  rapidly  de- 
graded due  to  the  effect  of  the  increasing  value  of  the  integral  term  in 
the  E  matrix. 
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4.  6   Estimate  of  Future  State 

As  noted  previously,  the  future  state  vector  of  a  continuous -thrust 
vehicle  is  not  a  deterministic  function  of  the  state  history  due  to  short 
period  random  variations  and,  perhaps,  long  term  degradation  of  engine 
performance.     If  engine  degradation,  presumably  in  the  form  of  loss  of 
specific  impulse,  does  not  become  a  factor,  the  space  navigator  would 
be  justified  in  assuming  a  zero  mean  for  engine  variations.     With  such 
an  assumption  a  filter,  biased  to  include  a  priori  information  of  engine 
statistics  and  a  priori  information  on  state  statistics,  may  be  construct- 
ed which  will  produce  an  optimum  prediction  of  future  state. 

If  long  term  variation  in  engine  performance  is  evident,  the  naviga- 
tor is  faced  with  the  problem  of  extrapolating  accumulated  engine  data 
in  order  to  make  a  prediction  with  acceptable  confidence       Certainly, 
sophisticated  analytic  techniques  exist  for  smoothing  and  extrapolating 
the  measured  engine  data.     However,  in  this  section  we  shall  only  be 
concerned  with  using  data,  regardless  of  the  manner  in  which  they  are 
processed,  to  predict  a  future  state.     For  example,  simple  graphical 
extrapolation  of  plotted  data  will  suffice 

Two  prediction  techniques  will  be  presented,   1)  a  very  simple  ex- 
trapolation technique  for  short  term  prediction  and  2)  an  optimum  filter 
technique  for  long  term  prediction.        The  first  method  assumes  the 
present  state  has  been  determined  by  the  methods  of  section  4   3  or  4.4 
It  uses  the  state  transition  equation  and  an  estimate  of  engine  perform- 
ance to  predict  state  in  the  near  future       The  estimate  is 

tn+1 

&s     .  -   =  T  x1       S£     +    A    "}.       f       AB  &£  dt  (4-25) 

—  n+1  n+ 1  n    —  n  n+ 1      J  — 

tn 

where  the  notation  t  =  t    ,  -  will  indicate  the  future.      This  method  may 

n+1  J 

be  characterized  as  an  extrapolation  of  filtered  data  and  may  be  quickly 
computed. 

The  second  method  uses  the  optimum  filter  of  section  4.  3  but  in 
this  case  a  representative  subvector  of  the  measurement  vector  is 
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6q  .  =  6r  .  +     /    D.  6a  dt  +    /    D.   6^  dt  (4-26) 

the  error  vector  associated  with  the  i      position  measurement  is 

ln  ln+l 

e.  =    e    .  +     f    D.e     dt  +     f  D.  £     dt  (4-27) 

—  l—  ri        J        l  —  a  ^       l  —  a  v  ' 
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where  6a  and  e      represent  the  extrapolation  of  engine  performance  and 

—  —  3, 

the  error  in  that  extrapolation,    respectively.  F 

The  optimum  biased  filter  of  section  4.  4  may  be  used  without  al- 
teration, however,  the  covariance  matrix  of  measurement  error  will  be 
more  difficult  to  compute. 

This  second  method  may  be  characterized  as  a  filtering  of  extrap- 
olated   data  and  should  be  used  for  long  range  prediction.     The  proof 

32 
by  Potter  and  Stern   'J  shows  this  second  method  to  be  optimum  and 

intuitively  it  appears  to  be  of  correct  form.     By  comparing  expressions 

for  the  covariance  matrix  of  uncertainty  in  state  at  time  t      .  for  the 

two  methods. 

U  =  <  u  uT>  (4-28) 

it  is  apparent  that  both  methods  result  in  the  same  error  ellipsoid  as 

t    ,  ,   approaches  t    .     For  both  methods: 
n+ 1     ^  n 

U  =  (NT  E"1  N)  (4-29) 

t    ,1 — «-t 
n+1  n 

thus  the  assertion  that  the  simple  technique  will  be  quite  adequate  for 
short  term  predictions  appears  justified. 
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CHAPTER  V 
TRAJECTORY  DETERMINATION 

5.  1   Summary  of  Chapter  V 

In  this  chapter  the  problem  of  finding  a  trajectory  which  satisfies 
the  end  condition.,  the  propulsion  restriction  and  which  minimizes  pro- 
pellant  consumption  is  discussed  and  illustrated  with  the  VSI  trajectory. 
The  problem  is  formulated  using  the  classical  methods  of  the  calculus 
of  variations.     A  new  interpretation  is  presented  for  the  results,  which 
are  shown  to  yield  the  same  control  program  as  the  Pontryagin  maxi- 
mum principle.     The  guidance  theory  of  Chapter  III  is  therefore  sug- 
gested as  a  trajectory  computation  scheme.      Finally,  qualitative  as- 
pects of  low-thrust  optimal  trajectories  are  discussed. 

5.  2    General  Remarks 


A  vast  amount  of  effort  has  been  directed  in  recent  years  toward 
the  study  of  optimization  problems.     Such  problems  belong  to  the  cal- 
culus of  variations  which  owes  its  early  development  to  such  men  as 
Lagrange,  Euler,  Hamilton  and  Gauss       The  introduction  of  high  speed 
computers  has  been  the  primary  impetus  in  bringing  renewed  interest, 
after  years  of  limited  application,  to  variational  techniques       Specifi- 
cally, much  recent  literature  treats  the  characteristics  of,  the  neces- 
sary conditions  for,   and  methods  of  computing  solutions  to  optimization 
problems.     Variational  techniques  are  used  almost  exclusively  as  the 
primary  mathematical  tool.     The  contributions  listed  as  references  in 
this  report  constitute  only  a  minor  fraction  of  the  published  works, 
These  efforts  have  resulted  in  a  large  body  of  theory  now  called  opti- 
mal control  theory  which  has  application  to  virtually  all  optimization 
problems  dealing  with  dynamical  systems. 

A  fundamental  precept  of  the  theory  is  that  along  optimal  trajec- 
tories, admissible  first  order  variations  in  an  unconstrained  control 
program  cannot  produce  a  first  order  effect  in  the  cost  function      For 
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unconstrained  controls,  admissible  control  variations  are  those  which 
do  not  cause  a  first  order  change  in  boundary  conditions.      This  princi- 
ple is  derived  directly  from  basic  equations  and  is  part  of  the  definition 
of  optimality  in  the  calculus  of  variations       A  mathematical  consequence 
of  this  principle  is  that  certain  arrays  of  coefficients,  when  evaluated 
along  the  optimal  trajectory,  will  have  a  zero  determinant.     When  work- 
ing with  the  practical  problem  of  computing  trajectories,  or  of  guidance 
around  the  optimal  trajectory,  the  inverse  of  a  matrix  defined  by  such 
arrays  usually  appears  in  the  equations.     Obviously  the  occurrence  of 
a  singularity  complicates  the  procedure  of  finding  an  optimal  control 

program  and  its  associated  optimal  trajectory.     It  has  led  to  investiga- 

39 
tion  of  second  variations  of  the  cost  and  the  trajectory      ,  to  various 

12 

schemes  for  finding  near-optimal  controls      ,  and  to  optimal  controls 

14 
which  only  approach  the  desired  boundary  conditions 

The  implication  of  these  investigations  is  disturbing  from  a  physi- 
cal viewpoint  since  they  imply  that  an  otherwise  well  behaved,  smoothly 
operating  system,  in  some  sense  becomes  uncontrollable  along  an 
optimal  reference  trajectory.     In  reality,  singularities  are  more  often 
mathematical  than  physical.     A  further  implication  is,  that  although  a 
control  can  be  found  which  approaches  the  optimal  control  to  within  a 
small  increment,  it  is  orders  of  magnitude  more  difficult  to  find  the 
exact  optimal  control. 

The  physical  world  does  not  usually  behave  in  such  an  unruly  man- 
ner, thus  the  answer  must  be:     1)  the  mathematics  have  an  interpreta- 
tion that  has  been  overlooked  or  2)  the  problems  may  be  approached 
from  a  different  viewpoint. 

Actually,  both  1)  and  2)  have  validity.     To  support  this  contention, 
the  problem  of  generating  an  optimal  reference  trajectory  for  use  in 
testing  the  guidance  theory  of  Chapter  III  is  considered  as  an  example 
To  be  sure,  a  conservative  field  such  as  the  gravitational  field  is  a 
well  behaved  space  in  which  to  work.     Undoubtedly,  problems  which 
deal  with  dissipative  forces  or  higher  order  nonlinearities  may  present 
difficulties  not  readily  resolved  by  the  method  of  this  chapter.     Hope- 
fully, however,  it  will  provide  an  approach  that  can  serve  as  a  starting 
point  for  more  difficult  problems. 
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5.  3    The  Calculus  of  Variations  Problem 

In  Chapter  II  it  is  shown  that  minimizing  the  acceleration  integral 
is  equivalent  to  minimizing  the  propellant  consumption  for  a  low-thrust 
vehicle.     Thus  the  mass  rate  equation  is  superfluous  for  VSI  trajecto^ 
riesanditis  necessary  only  to  work  with  the  equations  of  motion  and 
the  acceleration  integral.     For  the  VSI  vehicle  the  problem  can  be  for- 
mulated as:     Given 

r  (0)  =-r 
—  s    '      — o 


satisfy 


v  (0)  =  v 
—  v    '      — o 


£(*f)  =  If 


y(tf)  =  vf 


(5- 

-1) 

(5- 

■2) 

(5- 

-3) 

(5- 

4) 

subject  to 


r  =  y  (5-5) 

v  = r  +  a  (5-6) 

—  3    —      —  \        / 

r 

and  minimize  t^ 

J  =  —      /     a2  dt  (5-7) 

A  scalar  functional,,  F,  will  now  be  formed  using  the  well  known      :j      ; 
Lagrange  multiplier  technique. 

2 

F  =  X  T  (f  »  v)  +  X  T   (v  +  -%r    r  -  a)    +  —  (5-8) 

—  r     —      — '       —v     —         3  — '  n  ' 

r  I 

where  X      and  A      are  time  varying  Lagrange  multiplier  vectors  (or 
Euler  variables).     The  remaining  variables  have  the  same  meanings 
as  in  previous  chapters 

Since  the  bracketed  terms  are  zero,   clearly 

J  =     J     F  dt  (5-9) 

0 
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Integrating  by  parts  and  setting  the  first  variation  of  the  integral  equal 
to  zero,  one  obtains 


6    /     F  dt  =  0  = 


X       Sr  +  X    '   Sv 

—  r  —  v     — 


0 


f 


0 


•V 


+  (a      -  X     )  5a 
—  v 


—  v 


dt 


(5-10) 


Applying  the  fixed  end  conditions  and  setting  the  coefficients  of  the 
state  and  control  variations  equal  to  zero  gives  the  Euler  equations  and 
the  boundary  conditions. 
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(5-11) 
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i     t, 


(5-12) 
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J    0 


(5-13) 


A 


(5-14) 


Notice  that  the  Euler  variables,  (5-11),  satisfy  the  adjoint  relationship 
to  the  variational  equations  of  motion  and  that  the  boundary  conditions 
are  unknown.     Equations  (5-12)  and  (5-13)  are  written  to  emphasize 
the  unknown  boundary  conditions. 

Clearly,  if  the  correct  initial  value  of  the  six  component  Euler 
vector  were  known,  the  entire  system  of  state  and  Euler  equations 
could  be  integrated  simultaneously  from  t  =  0  to  t  -  tf.     The  unique 
reference  trajectory  and  its  associated  optimal  control  program  would 
then  be  known.     A  procedure  for  finding  the  initial  value  is  discussed 
in  later  sections 
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A  very  annoying  fact  about  the  preceding  approach  to  finding  ex- 
tremals in  the  calculus  of  variations  is  the  necessity  of  searching  for  a 
formulation  of  the  problem  which  gives  a  meaningful  result      Minimizing 
the  cost  function  J  is  known  to  be  equivalent  to  minimizing  propellant 
consumption,  thus  one  differential  equation  of  state  was  eliminated,  in 
particular,  the  differential  equation  describing  the  rate  of  change  of  the 
optimized  quantity,  mass.     A  simple  answer  resulted.     It  should  not  be 
necessary  to  search  for  an  equivalent  formulation  for  a  problem  if  the 
set  of  differential  equations  describing  the  system  is  complete,  linearly 
independent  and  reasonably  well  behaved.     Two  different  Mayer  formula- 
tions of  the  VSI  problem  help  illustrate  this  argument. 

For  the  first  problem,  minimize  propellant  consumption  instead  of 
the  acceleration  integral.     Equations  (5-T)  through  (5-6)  hold,  but  the 


mass  rate  equation  must  be  reintroduced.     Thus  it  is  required  that 

m=-*?™i  (5-15) 

2p 

and  that 
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m   dt  =  m 
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m, 


(5-16) 


be  a  minimum  where  m     is  propellant  mass.      For  fixed  initial  mass 
the  first  variation  of  m     equals  the  variation  of  -  mf ,      Forming  the 
functional  F 


T     •  T     •        u 

F  =  A       (r  -  v)  +  A       (v  +  -4  r 
—  r     —      — '       —  v     —         3   — 
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2      2 
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(5-17) 

then  integrating  by  parts  and  setting  the  first  variation  of  the  integral 
equal  to  zero  as  before,  yields 
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In  this  case  the  Euler  variables  for  position  and  velocity  again  satisfy 
equations  (5-11)  through  (5-13).     In  addition 


(A       -  1)     -     ?  (5-19) 

0 


A        =   1  (5-20) 

m, 


f 


2 

A       -  A        ^L     --  0  (5-21) 

mm  v  ' 


A 
a  =   ~-  (5-22) 

A       iB_ 
m 


Equation  (5-22)  appears  to  differ  from  (5-14)  by  a  scalar  function  of 
time.     Since  the  control  program  is  unique  (5-22)  must  reduce  to 
(5-14).     By  manipulating  (5-15)  and  (5-22)  it  is  possible  to  obtain 


A 


=   -  2  —  (5-23) 
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which  has  the  solution 


%        .  constant  ^-24" 

m  ™  2 

m 


Applying  the  boundary  condition  (5-20) 

2 


Substituting  into  (5-22) 


mf 

A       -   — —  (5-25) 

m  2 

m 


a  =  4"^v  (5-26) 

mf 
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Equation  (5-26)  differs  from  (5-14)  only  by  a  constant  which  is  easily 

absorbed  in  the  unknown  boundary  conditions  on  A    .     The  Euler  variable 

J  —  v 

A      appears  to  have  been  superfluous. 

Consider  a  slightly  different  formulation  which  is  just  as  valid.    Use 
thrust  as  the  control  then 

r  =  v  (5-5) 

u  f 

v  =  _  JL    r  +  —  (5-27) 

—  3    —      m  v  / 

r 

f2 
m  =  -  —  (5-28) 

2p 
minimize  m     =      J      -  m  dt  (5-2  9) 

Proceeding  as  before,,  the  functional  F  and  its  first  variation  are 

f  2 

F  =  AT(r-v)  +  AT(v  +  -^rr-  —  )+A       (m  +  —  )  -  m 

-r   V_      _>      _v   *_      r3   -      m>         m  2p 


to  t  p  t  £• 


(5-30) 
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(5-31) 


m  p 
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The  blank  brackets  contain  the  same  terms  as  (5-18)  and  are  used  to 
show  the  differences  between  the  two  formulations.     In  this  case 

T  -f 

\n  -  iv  "I  (5"32) 

m 

A 

(5-33) 


A      ™ 

m   p 


63 


T 
Using  (5-33)  to  eliminate  X      from  (5-32)  and  using  the  state  equation 

(5-28),  it  is  found  that  again 

A       =  -  2  X     —  (5-23) 


m  m 


m 


Solution  (5-25)  holds  thus 

f  =-£-*    m  A  (5-34) 

—  2  —v  v 

mf 

which  is  the  same  as  (5-26)  since  a  =  f/m. 

The  preceding  VSI  examples  result  in  the  same  control  and  in  none 
of  them  is  it  necessary  to  use  the  variable  X     .     Before  interpreting 
what  this  means  one  additional  example  is  presented  which  involves 
placing  a  constraint  on  the  thrust.     Assume  that 

f2  <    f2  (5-35) 

—      o  \  / 

17 
where  f     is  the  thrust  limit.     Following  Breakwell      ,  the  constraint 
o  & 

may  be  handled  by  adding  to  the  functional  (5-30)  the  term 

7     (f2  -  f2)  =  0  (5-36) 

where  y    is  selected  to  satisfy  (5-36).     Thus  for  this  case 

F  =  [(5-30)J+/  (f2  -  f2)  (5-37) 
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where  the  blank  brackets  are  the  same  as  (5-31)  and  not  needed  in  the 
argument. 

Consider  the  coefficients  of  6m  and  6f  in  (5-38) 

T     1 
A       -  A  l *-    =  0  (5-39) 
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f  = 
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—  v 
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(5-40) 
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Eliminating  A       in  (5-39)  by  using  (5-40) 
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11 
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(5-41) 


When  f  =/  f  ,  y     must  be  zero  to  satisfy  (5-36).     Therefore  the 

solutions  for  A       and  f  in  such  unconstrained  regions  can  differ  from 
m  —  to 

(5-25)  and  (5-34)  at  most  by  a  constant      When  f  =  f     the  only  admissi- 
ble solution  for  (5-40)  is 
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Therefore 
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With  this  result  it  is  possible  to  rewrite  (5-40)  with  a  new  variable  y 
such  that 


f  =  my  A 
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(5-44) 


where 
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The  Euler  variable  A       and  its  differential  equation  have  again  been 
eliminated  from  the  problem. 

In  the  preceding  examples  the  manipulations  required  to  remove 
A       are  not  particularly  difficult  but  by  no  means  are  they  obvious.     The 
meaning  of  a  superfluous  Euler  variable  immediately  comes  into  ques- 
tion and  it  seems  reasonable  to  inquire  if  there  is  general  significance 
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in  this  phenomenon  or  if  the  preceding  examples  are  only  isolated 
special  cases.     Conceptually  there  is  nothing  about  the  variable  A 
which  is  unique  except  for  its  association  with  the  cost  variable  m  in 
a  Mayer  form  problem      In  reference  17  Breakwell  interprets  the 
Euler  variables  as  the  sensitivity  of  the  absolute  minimum  of  the  cost 
function  (i.  e.    the  unconstrained  minimum)  to  changes  in  the  state  vari- 
ables.    Therefore  for  the  problem  of  this  thesis 

T           3  m''              3  m! 
A    (t)  =  2_    = L  (5-46) 

8^t  8*t 

where  m"''  denotes  the  minimum  value  of  m     obtainable  from  the  initial 

P  P  * 

state  when  no  control  constraints  are  used;  similarly  for  -  mf . 

In  a  typical  two  point  boundary  value  problem  with  physical  rates 
of  the  form  s  =  g  (s„a),  the  usual  objective  is  to  find  a  control  solution 
which  satisfies  the  end  point.     If  any  of  the  state  variables  have  free 
end  conditions  the  solution  is  not  unique.     Optimization  criterion  are 
then  used  to  assure  uniqueness.     From  Breakwell's  interpretation, 
when  the  Euler  variables  are  formed  into    i  vector  they  describe  the 
direction  in  state  space  of  maximum  sensitivity  of  the  cost,  i.  e.    the 
gradient.     However  for  state  variables  involved  in  the  cost  function 
such  information  is  available  directly  from  the  governing  differential 
equations.     In  answering  the  question:   "How  do  state  variations  affect 
the  cost?  ",  the  Euler  variables  provide  a  coupling  between  state  vari- 
ables and  the  cost  function  which  for  some  variables  is  not  evident 
from  the  system  equations.     But  for  state  variables  involved  in  the 
cost  this  coupling  is  essentially  a  priori  information  and  is  in  some 
sense  superfluous. 

The  classical  approach  in  the  calculus  of  variations  assigns  an 
Euler  variable  for  each  state  variable.     If  a  variable  is  superfluous  it 
should  drop  out  of  the  formulation.     But  finding  a  way  to  assure  that  it 
does  may  be  exceedingly  difficult,  even  when  it  is  recognized  that  the 
variable  is  superfluous.      As  is  subsequently  shown,  the  Euler  variables 
correspond  to  the  costate  variables  in  Pontryagin  formulations.     If  the 
variable  is  superfluous  in  the  calculus  of  variations  it  is  also  superfluous 
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using  the  Hamiltonian  approach.     This  fact  was  tacitly  acknowledged  in 
Chapter  III  when  the  state  variation,  5mf,was  dropped  from  considera- 
tion and  the  last  row  of  the  adjoint  set  was  deleted      The  reason  given 
in  Chapter  III  for  deleting  one  row  of  adjoint  functions  was  different 
from  here,  but  the  result  is  the  same. 

The  general  significance  of  a  superfluous  Euler  variable  is  indi- 
cated by  comparing  the  classical  approach  to  the  newer  state  space 
approach  using  the  Pontryagin  maximum  principle.      The  classical 
formulation  suffers  from  several  difficulties,,  two  of  which  are:     1) 
The  number  of  variables  and  equations  are  usually  quite  large  and, 
perhaps  due  to  tradition,  are  often  treated  individually  as  scalar s      As 
a  result  algebraic  detail  often  obscures  important  ideas  in  classical 
formulations.     The  more  compact  vector  and  matrix  notation  are  only 
now  becoming  widespread  in  the  literature.     2)  The  theory  offers  no 
suggestion  of  how  to  solve  for  the  correct  boundary  conditions  on  the 
Euler  variables, 

However,  an  important  attribute  of  the  more  detailed  notation  is 
indicated  by  the  preceding  examples       By  sufficient  manipulation  an 
unneeded  variable  was  eliminated      It  is  more  than  coincidence  that 
the  variable  is  precisely  the  one  that  creates  problems  in  the  state 
space  approach 

The  state  space  approach  using  vector  and  matrix  notation  as  in 
Chapter  III  is  popular  because  of  its  compactness       Further,  when 
coupled  with  the  Pontryagin  maximum  principle  it  may  be  used  to 
solve  optimization  problems  including  cases  of  linear  control  (the  CSI 
problem)  which  can  not  be  completely  solved  in  the  calculus  of  varia- 
tions.    The  approach  also  provides  a  method  of  determining  the  bound- 
ary values  for  the  costate  variables.     Finally,  the  optimality  criterion 
is  more  useful.     This  last  statement  is  discussed  fully  in  Appendix  D. 

In  the  state  space  approach  the  system  differential  equations  are 
written  in  vector  form,  an  appropriate  cost  function  is  chosen  and  then 
methods  are  sought  to  find  the  desired  solution.     From  the  references 
previously  mentioned  it  is  apparent  that  the  solutions  often  encounter 
a  singular  matrix      In  the  problem  of  this  thesis  the  singular  matrix 
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is  denoted  as  M     and  is  the  seven  by  seven  matrix  corresponding  to 
M  in  Chapter  III.     These  matrices  are  discussed  in  detail  in  Appendix  C. 
It  is  suggested  that  such  singular  matrices  may  be  the  result  of  attempt- 
ing to  optimize  a  function  of  state  and  at  the  same  time  retain  a  costate 
variable  (or  Euler  variable)  which  is  superfluous.     In  the  problem  of 
this  thesis  such  was  the  case.     It  is  further  suggested  that  difficulties 
sometimes  encountered  in  numerically  finding  terminal  values  for  the 
Euler  variables  may  also  be  due  to  the  presence  of  a  superfluous  Euler 
variable  which  was  not  recognized  as  such. 

No  attempt  will  be  made  in  this  report  to  derive  the  general  rule 
which  will  show  when  an  Euler  variable  is  superfluous  and  thus  may  be 
eliminated  a  priori.     It  may  be  the  case  that  those  matrix  elements 
which  cause  nonphysical    singularities  to  occur  in  state  space  formula- 
tions can  always  be  removed  without  changing  the  problem. 

In  the  remainder  of  this  chapter  the  relationship  of  the  calculus  of 
variations  solution  to  the  state  space  solution  is  shown  for  the  low - 
thrust  problem  only.     The  method  of  deleting  A.      is  illustrated. 

5.  4   Removal  of  a  Superfluous  Euler  Variable 

Although  an  Euler  variable  may  be  superfluous,  a  general  method 
of  removing  one  is  not  obvious  from  the  derivation  in  section  5.  3.     Con- 
sider, however,  the  following  derivation  of  the  problem  in  section  5.  3 
which  combines  the  compact  notation  with  the  classical  approach.     It 
may  be  indicative  of  a  general  method. 

Assume  the  system  is  described  by  physical  rates  which  have  the 
form 

L  =  I  (£*  a)  (5-47) 

where  s  is  the  state  vector,  a  the  control  vector  and  g  is  a  vector  func- 
tion of  state  and  control. 

Assume  that  the  cost  variable  is  some  scalar  function  of  the  physi- 
cal rates  or  some  nonlinear  function  of  the  control. 

S  =  h  (s)  (5-48) 
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or 

S  =  h  (a)  (5-49) 

where  S  is  the  cost  integrand 

Further,  allow  the  physical  rates  to  be  subjected  to  constraints. 
Define  the  constrained  variables  as  a  vector  function  of  the  physical 
rates  or  of  the  control.     That  is 

Z  .  =  Z  [£  (s,  a)]  (5-50) 

A  general  functional,  equivalent  to  those  used  in  section  5.  3  is 

F  =  S  +  ATg  +  y'Ty  (5-51) 

where  A  is  the  vector  of  Euler  variables  and  where  y    is  determined 

in  such  a  fashion  that  the  constraint  is  satisfied.     It  is  a  vector  equiva- 

i 
lent  of  the  y     in  section  5   3 

Application  of  the  usual  variational  technique  produces  the  expanded 
set  of  state  equations,  Euler  equations  and  control  equations, 

s  =  g  (s,  a)  (5-47) 

-mX-(°k    +XT    »I„+yTd-l)T  (5-52) 


^3s  3s  3s 

0=   m    +AT    !i    +Y'TdL    f  (5-53) 

^3a  3a  3a    / 

For  the  unrestricted  case  such  that  y  =  O  and  using  the  familiar 
definitionof  the  matrices  A  and  B,  that  is 

£i  =    A  (5-54) 

3s 

^1  =  B  (5-55) 

3a 


then  from  (5-52)  one  obtains 


T  T 

A      =  -  A      A  (5-56) 
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A  solution  to  the  set  (5-53)  to  within  an  unknown  scale  factor,  is 

a  =  BT  X  (5-57) 

That  (5-57)  is  a  solution  may  be  verified  by  expansion.     The  variational 
equations  of  state  and  their  adjoint  set  are 

6s    =  A  6s  +  B  6a 


(5-58) 


A=  -Aa 


(5-59) 


Since  the  Euler  equations  (5-56)  satisfy  the  adjoint  relation  as  does 
(5-59)  a  solution  for  (5-56)  is  obtained  by  choosing 


Af  =  i 


Then 


Ai% 


(5-60) 


(5-61) 


Dropping  the  subscript  t  from  (5-61)  and  substituting  into  (5-57)  yields 


a 


BT    ATA 


(5-62) 


By  partitioning  TV   and  X  f  it  is  apparent  that  deleting  the  seventh  row 


f 

of  TV  3     A.  n  ,  and  deleting  the  unknown  boundary  value  X        eliminates 

—  7  b  J  mf 

the  superfluous  Euler  variable  from  consideration.     That  is 

f  A 


.  -  b*  [A*T  j   A7] 


■v 


(5-63) 


m 


f 


In  order  to  find  the  optimal  trajectory  it  is  only  necessary  to  find 
the  final  value  of  the  six  Euler  variables.      (Or,  since  the  adjoint  and 
fundamental  solutions,  discussed  in  Chapter  III,  allow  transformation 
at  will  between  one  terminus  and  the  other;  the  initial  values  of  the 
Euler  variables  may  be  used.  ) 

By  letting 

(5-64) 
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and  deleting  -A_     and  A        from  (5-63)    ^h0  expression  yie 


Ids 


a  =  BT  A    '    v  (5-65) 

which  is  identical  with  the  expression  for  the  optimal  control  program 

for  guidance  derived  in  Chapter  III.     Thus  both  the  Pontryagm  principle 

and  the  Euler  variables  give  the  same  result  when  A      is  deleted.      It  is 

&  m 

shown  in  Appendix  C  that  along  the  optimal  trajectory  Af  is  a  null  vec- 

T     A  T  ~ 1 

tcr  of  B     J\    ,  thus  equation  (5-62)  displays  the  well  known  and  vei  v 

troublesome  singularity  when  evaluated  along  the  optimal  trajectory. 

This  is  the  consequence  of  optimality  which  has  instigated  the  search 

for  new  methods, including  this  one.     It  is  further  shown  that  v_  is  not  a 

null  vector  of  B     A.        thus  (5  -65)  can  be  use   ' 

| 

For  the  re^+r:cted  VSI  case  the  preceding  arguments  hold  wi+h  y 
entering  as  in  section  5,  3.     The  restricted  VS    i     se  is  ccmp>  I 
determinate  in  the  calculus  i  *  vai  Lations  using  Breakwell  s  approach 
and  deleting  A     „     The  CSI  case  may  also  b<=  treated  in  an  analogc 
fashion  except  that  the  coast  phase  can  n<  '   '        niquely  determined. 

This  is  due  to  the  so-called  degeneracy  of  the  calculus  of  va]      fci< 

7 
for  linear  optimums    . 

5.  5   Solution  by  Direct  Integration 

There  are  two  general  approaches  for  finding  an  optimal  trajectory. 
Descriptions  of  the  first  me*hcd  are  usually  prefaced  by:  "If  any  non- 
optimal  solution  can  be  found  which  satisfies  the  boundary  conditions., 
then  the  solution  can  be  moved  in  th<    direction  of  the  optimal,  etc,  " 
This  approach  is  usually  called  a  gradient  meth<  ■   :  -  re1  used  m 

this  report  because  nonoptimal  solutions  cannot  satisfy  the  Euler 
equations,  thus  it  would  be  difficult  to  relate  any  nonoptima]  c<  ntr< 
to  the  form  of  equation  (5-65). 

The  second  general  approach  is  often  described  as  "solving  \he 
wrong  problem  in  an  optimal  manner.  '     This  is  the  approach  used  in 
Chapter  VI      An  initial  guess  is  made  for  the  initial  Euler  variables 
Vj,  then  the  state  equations,  adjoint  set  and  Euler  variables  are  integrat- 
ed to  the  final  time.     The  resulting  final  values  are  compared  with  the 
desired  final  values  and  a  new  es'   m    te  of  v  is  computed.     The  process 
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is  repeated  until  the  boundary  values  are  satisfied.     The  fact  that  the 
gravitational  field  is  conservative  and  well  behaved  aids  materially  in 
improving  the  speed  of  convergence. 

5.  6   Properties  of  Optimal  Trajectories 

It  is  possible  to  gain  considerable  insight  into  the  properties  of 
propellant-optimal  trajectories  without  resorting  to  machine  computa- 
tion.    Such  insight  aids  materially  in  working  with  the  mathematics  of 
the  transfer  process  and  often  leads  to  new  methods  of  approach. 

For  this  purpose,  consider  the  problem  of  transferring  between  two 
planets  with  minimum  propellant  expenditure  as  one  of  changing  energy 
and  angular  momentum  in  the  most  efficient  manner. 

From  the  energy  integral  of  orbital  mechanics  one  may  write  the 

energy  per  unit  mass  of  the  vehicle  as 

2 
e  =  —  -    -^  (5-67) 

2  r 


and  the  energy  time  derivative  as 


v  v     + 


-^2     r  (5~68) 

r 


Equation  (5-68)  may  be  written  in  the  vector  notation  and  reduced  with 
help  of  the  equations  of  motion  to 

e    =  yT  a  (5-69) 

Equation  (5-69)  represents  the  rate  of  energy  change  imparted  to  the 

vehicle.      The  rate  of  energy  expenditure  by  the  propulsion  system  is 

the  power 

•      2  •      T 

m  c  m  c      c  ,r    7m 

p  = —    =    -  — = — =-  (5-70) 


where  rh  is  considered  as  the  mass  flow  rate  per  unit  mass  and  £  has 
the  opposite  sense  of  a. 

Similarly  a  functional  expression  for  angular  momentum  is 

h  =  r  X  v  (5-71) 
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h=r.Xy+rXy  (5-72) 

where  the  notation  signifies  the  vector  product  of  r  and  v.     Equation 
(5-72)  is  reduced  with  help  of  the  equations  of  motions  to 

h  =  r  X  a  (5-73) 

Using  equations  (5-69)  and  (5-73)  the  problem  is  formulated  as: 

tf 

Minimize  S  =     /     -  m  dt  (5-74) 

0 


subject  to  the  constraints 


tf    T 
Ae  -     /  v     a  dt  (5-75) 

0 

*f 

Ah  -     /  i    <  a  dt  (5-76) 

0 

plus  velocity  and  position  constraints  and  where   Ae  and  Ah  are  known 
quantities. 

The  above  formulation  is  not  computationally  desirable:  however 
it  provides  a  basis  for  seme  significant  qualitative  arguments. 

It  is  apparent  from  equation  (5-69)  that  with  a  fixed  in  magnitude 
the  vehicle  changes  energy  most  rapidly,  thus  most  efficiently,  when 
v  is  large  and  y  and  a  are  colinear,      From  the  energy  equation  (5-67) 
one  observes  that  in  a  central  force  field  the  velocity  magnitude  goes 
inversely  with  the  square  root  of  the  radius.      From  equation  (5-73)    it 
is  apparent  that  h  changes  most  rapidly  when  r  is  large  and  orthogonal 
to  a. 

From  the  preceding  observations  one  may  deduce  that  the  trajec- 
tory which  minimizes  the  propellant  expenditure  will  tend  to  keep  the 
exhaust  velocity  vector  aligned  with  the  large  velocity  vector  when  the 
vehicle  is  deep  in  a  gravitational  field  so  that  energy  is  changed  most 
efficiently.     It  will  also  tend  to  rotate  the  angular  momentum  vector 
when  the  vehicle  is  far  out  in  the  gravitational  field  so  that  the  vector 
r  has  large  magnitude.     Where  these  requirements  are  incompatable 
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with  boundary  conditions,  the  optimum  solution  should  resolve  the 
conflict  in  favor  of  propellant  conservation. 

The  trajectories  which  have  been  generated  numerically  in  this 
thesis  all  appear  to  satisfy  these  precepts  in  so  far  as  boundary  con- 
ditions permit.     The  increase  in  specific  impulse  which  is  character- 
istic of  variable  thrust  rockets  as  they  traverse  the  center  portion  of 
the  heliocentric  phase  is  due  to  the  acceleration  vector  becoming  orthog- 
onal    to    the    velocity  vector.     Likewise  the  coast  phase  for  a  thrust- 
limited  rocket  occurs  when  the  acceleration  vector    rotates   through 
the  orthogonal  orientation. 

With  the  qualitative  arguments  of  this  section,  it  is  possible  to 
sketch  a  reasonable  approximation  to  a  low-thrust  optimum  trajectory 
without  resorting  to  machine  computation. 
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CHAPTER  VI 
COMPUTATIONAL  PROCEDURE  FOR  TRAJECTORIES 

6.  1   Summary  of  Chapter  VI 

Details  of  a  procedure  for  using  the  guidance  equation  as  the  basis 
for  trajectory  computations  are  set  forth  in  this  chapter.     Specifically,, 
the  differential  equations  to  be  integrated  and  the  procedure  for  cor- 
recting initial  conditions  on  the  Euler  equations  are  presented. 

6   2    General  Remarks 


In  Chapter  V  the  fundamental  problem  of  selecting  initial  (or  final) 
conditions  on  the  Euler  variables  is  introduced.     Since  the  acceleration 
program  is  an  explicit  function  of  the  Euler  variables  it  is  necessary  to 
make  a  first  estimate  of  the  Euler  initial  conditions  in  order  to  start 
the  iteration      This  estimate  will  usually  be  a  gross  error  under  the 
best  of  circumstances,  thus  it  is  desirable  to  choose  a  very  simple 
estimate  such  as  the  null  vector  or  a  unit  vector.     Because  of  this,    the 
procedure  we  seek  must  be  strongly  convergent  and  independent  of  the 
error  in  the  first  iteration  attempt. 

In  any  iterative  procedure  the  speed  of  convergence  is  directly 
dependent  upon  the  validity  of  assuming  that  certain  quantities  do  not 
change  from  one  iteration  to  the  next.     In  the  preceding  chapters  it  was 
explicitly  assumed  that  the  miss  vector  at  the  terminal  point  was  de- 
pendent only  upon  a  change  in  the  initial  Euler  vector  and  that  A.  and 
B  were  invariant  for  neighboring  trajectories.     For  small  perturbations 
around  a  reference  trajectory  such  an  assumption  is  valid..     However 
for  the  large  perturbations  expected  in  trajectory  computation  A  and 
B  are  not  invariant.     The  change  of  these  quantities  from  one  iteration 
to  the  next  and  the  change  of  other  second  order  quantities  will  deter- 
mine the  speed  of  convergence  of  the  procedure 
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6.  3    The  Correction  Procedure 

From  Chapter  III  one  may  write  the  optimal  acceleration  program 
as  a  scalar  times  a  matrix  function  of  the  adjoint  set  and  the  initiali- 
zation vector. 

o         o  t-,  T    a  *  T     o  ,e    a  \ 

a     =  y     B      /V        v  (6-1) 

where  y  is  the  scalar  switching  function  derived  from  the  limiting  con- 
dition on  acceleration,  B  is  the  seven  by  three  matrix  of  partial  deriva- 
tives, 3g/3a,    A.*  is  the  six  by  seven  reduced  adjoint  set  and  v_  is  the 
six  by  one  initialization  vector. 

In  order  to  simplify  the  notation  in  certain  equations  to  follow, 
Euler  variables  will  be  used  interchangeably  with  the  matrix  notation. 
Recall  from  Chapter  V  that 


and 


A      =  BT  A*T  v  (6-2) 

—  v  — 


A      =   -  A  (6-3) 

—  v  —  r  ' 


The  objective  of  the  iterative  search  for  the  optimal  acceleration 
program  is  to  find  the  vector  v_  which  causes  the  six  component  miss 
vector,  J_,  to  approach  the  null  vector  to  some  desired  accuracy.     Thus 

v      =  V       1  +  Av  (6-4) 

—  n      — n-1  — 


and  we  desire  that 


v      =   v°  (6-5) 

—  n      —  v         ' 


for  n  as  small  as  possible,  where  n  is  the  iteration  number. 

The  miss  vector,  |  ,  may  be  written  as  the  difference  in  final  state 
which  results  when  a  nonoptimal  acceleration  program,  a,  is  used  in 


place  of  a   .     Define 


77  =     J        A*B  a  dt  =  Mv  (6-6) 

0 


and  f 


|  =     /    A*  B  (a°-a)  dt  (6-8) 

0 


76 


then  using  (6-1)  and  (6-8) 


|  =  /  A  *  B  BT  A  *T  *(yv)  dt 

0 


(6-9) 


Equation  (6-9)  implicitly  assumes  that  A  and  B  are  invariant  functions 
of  time  and  thus  constant  from  one  iteration  to  the  next       The  assump- 
tion permits  adequate  convergence  and  is  used  throughout. 

In  order  to  provide  for  the  discontinuities  which  occur  in  certain 
variables  at  the  beginning  and  end  of  restricted  thrust  regions,  it  is 
convenient  to  write  equation  (6-9)  as 


J-l 


t,- 


i+1 


i=0   du     U 


y  v  dt    &v 


(6-10) 


where  j  is  the  number  of  distinct  phases  of  the  trajectory      Expanding 
(6-10)  according  to  Leibniz's  rule  yields 


-i-  S 

i=0 


r-  t 


i+1 


L*l 


/  (A  b  BTA*T  y  +  A:  B  BTA*Ty  ^-)  dt 

'  dp 


+  A  *  b  btA*  T  y  v 


*i+i   9+ 


"i 


A  v 


(6-11) 


(k  =  i,  l+l) 

Since  y  is  the  discontinuous  switching  parameter,  the  time  at  which  a 
discontinuity  occurs    t,     and  its  derivative  with  respect  to  v  may  be  eval 
uated  from  the  derive    Ion  in  Chapter  IV.     Rewriting  equation  (3-76)  we 
obtain 

— -  (6-12) 

dq 


uk 
dp 


dp 


at 


By  comparing  (3-52)  and  (6-2)  observe  that 


(6-13) 


77 


Using  (6-2)  and  (6-3)  and  taking  the  derivative,  yields 


at, 


dv 


A  T  BT  A*T 

—  V 

T 
A        A 

—  v     — r 


(6-14) 


The  term  dy  j  dv_  may  be  evaluated  from  the  definitions  in  Chapter  III, 
which  are  written  here  in  more  convenient  form. 

a 


y   -  mm 


1; 


o  1 


m        A 

—  v 


(VSI) 


(6-15) 


2a 


0     if 


o  1 


m 


>    1 


A 


y  =  < 


V 


>  (CSI) 


o        1 
m 


2a  - 

if  — ^    .    *    .    <    1 


A 


m 


■v 


A 


■v 


(6-16) 


Therefore  by  differentiating  and  applying  equation  (6-2) 


JBTA *t 
a       A       B    _/Y 
o     —v 


®2i     =< 

dv 


m 


0 


A 


v 


if  y  1  °,  i 


if  y  =  0,    1 


(6-17) 
From  equation  (6-6)  let  us  make  the  following  definitions  for  77   and  M. 

(6-19) 


A*  B  BTA*T 


77    =  y  J\.'f 

77   =  M  v_ 
Using  (6-11)  and  (6-14)  through  (6-20)  one  obtains 


J  (  ■       y  v  v 
0     v 


a    v  2 


dt  Av 


V  V 


T 


(£) 


I     T 
7  A      A 

'    —v  — r 


where  we  let 


y 


,      fo      7=  i,  ol 

/1,  0  j 


(6-20) 


'cut  off  -  e 
Av_ 

cut  on  +  e 

(6-21) 


7  = 

7 


(6-22) 
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The  final  term  in  equation  (6-21)  is  evaluated  at   points   such  that  y  is 
not  zero,  i.  e„   e    before  cut-off  and  e    after  cut -on. 


Define 


M*  =      f    f M 
0       \ 


then 


y  77  7] 


(a     \  t 


Au   -- 


dt  - 


.     T 

7]  7] 


I     T 

y    A       A 
•     _  v    _r 


off 


on 


(6-23) 


M- 


: 


and  from  (6-4),  (6-5)  and  (6-6) 

M*  "  l  | 
Thus  the  correction  procedure  to  be  employed  is 


o  -1 

v        M     ?) 


— n      — n-1 


M       L  k 


(6=24) 


(6-25) 


16-26) 


6„  4   Quantities  to  be  Computed 

In  order  to  obtain  the  quantities  M,;    and  I  of  equation  (6    26)  it  is 
necessary  to  integrate  the  state  variables,  the  adjoint  set,  the  vec4<  l 
fl  and  the  matrix  M*„     In  addition  1+  is  convenient  to  integrate  the  Euler 
variables  rather  than  compute  them  from  other  quantities. 

From  the  discussion  in  Chapter  V  it  may  be  noted  that  the  adjoint 
equations  are  to  be  integrated  backward  in  time  frcm  Af  =  I.     However, 
i+  is  easily  shown  that  integrating  backwards  in  time  is  an  unnecessary 
complication.     In  Chapter  III  the  variational  equations  of  state  and  their 
adjoint  set  were  solved  together  to  produce 


Af  %f =  At  6-t  +  -f  ^-B  6-  dt 


»6-27) 


and  A    was  arbitrarily  chosen  equal  to  the  identity  in  the  guidance 
problem.     For  computing  the  entire  trajectory  from  t  -  0    to  t  =  tf  one 
may  write  an  equivalent  equation 


A*  *  |  =  A  *  6S      +     /  A"   B  B  -    A    j  A  (y  V)  dt 
f  o     -o      0J 


(6-28) 
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where  Af     has  both  the  last  row  and  last  column  of  A  deleted  in  order 
to  be  compatable  with  ^  and  where  6s      =  O  since  the  problem  starts 
from  a  given  state.     Since  the  initial  condition  on  TV  is  still  arbitrary 
it  is  convenient  to  choose  _A_     =  I  so  that  backward  integration  is  not 
necessary.     Thus  in  place  of  (6-24)  use 

-A **  i  =  M*   Au  (6-29) 

A     -  I  (6-30) 

o  v  ' 

An  interpretation  of  this  procedure  is  as  follows.     An  arbitrary  accel- 
eration program  is  selected  (i.  e.    estimate  v  )  ;  then  the  state  and  ad- 
joint equations  are  integrated  forward  resulting  in  a  miss  in  position 
and  velocity,   -  £  .     Since  the  acceleration  program  satisfies  the  condi- 
tions for  an  optimum  to  whatever  terminal  point  it  reached,  namely 

(r 

the  propellant  required  to  reach  that  point  is  minimum. 

Thus  we  are  assured  that  when  |  -*-  O,  iru— ^-  iru.  .  also.     It  is  un- 

—        —        f  f(max) 

important,  therefore,  that  the  "miss"  in  mf  is  unknown,  since  if  the 
procedure  converges  on  the  physical  boundary  conditions  it  does  so  via 
a  propellant-optimal  trajectory. 

The  term  -  Af      £,  merely  reflects  the  physical  miss  back  to  the 
initial  state  and  compensates  for  the  procedure  of  integrating  the  ad- 
joint set  forward  in  time  instead  of  backward  in  time.     In  the  machine 
procedure  it  is  convenient  to  use  equation  (6-31)  in  lieu  of  (6-26). 

v  '   =  v       1   -  M*~y\_  t*|  (6-31) 

—  n      — n-1  f    —  \  i 

The  full  set  of  equations  that  must  be  integrated  are  listed  with  the 
appropriate  initial  conditions. 

s   =  g  (s,  a)  (6-32) 

s  (0)    =  sq  (6-33) 

A=  -  A  A  (6-34) 
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A     =  i 

o 


A.  =  -  Ax  A 


A  (0)  =  v 


A  *  BBTA 


ta*t   y  nn 


T 


M     =  y 


M     (0)  =  O 

U  =  M  K 
where  M  is  the  first  term  of  M 

M  (0)    =  O 


a    \  2 


m 


(6- 

■35) 

(6- 

-36) 

(6- 

■37) 

(6 

■38) 

V  (0)    =  o 


(6-39) 
(6-40) 

(6-41) 
(6-42) 


An  information  flow  chart  and  the  FORTRAN  program  used  to  mecha- 
nize the  precedingequations  are  included  in  Appendix  G. 

The  procedure  worked  satisfactorily,  converging  to  very  small 
values  for  j;  with  surprisingly  few  iterations.  The  initial  estimates 
for  v  that  were  used  are 


v  =  O 


(VSI  trajectories)         (6-43) 


°3 


v  =  { : 


0 


)        (CSI  trajectories)         (6-44) 
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CHAPTER  VII 
A  NUMERICAL  EXAMPLE 

7.  1   Summary  of  Chapter  VII 

In  this  chapter  a  sample  heliocentric  transfer  from  Earth  to  Mars 
is  described.     The  selected  mission  was  programmed  for  machine  com- 
putation to  test  the  theories  developed  in  the  thesis.     Numerical  results 
from  the  test  are  presented  and  discussed. 

7.  2    General  Remarks 


The  linear  guidance  theory  of  Chapter  III  assumes  that  an  optimal 
reference  trajectory  is  known.     In  order  to  test  the  guidance  theory, 
therefore,  it  is  necessary  to  compute  an  optimal  reference  trajectory. 
In  Chapter  V  the  discussion  of  a  superfluous  Euler  equation  shows  that 
the  solution  of  the  linear  guidance  equation     satisfies  the  necessary 
conditions  for  an  optimum.     Thus  in  theory,  iterative  solution  of  the 
guidance  equation  will  generate  the  desired  optimal  reference.     If  the 
technique  converges  the  linear  guidance  theory  is  validated.     If  the 
procedure  converges  rapidly  then  the  theory  also  provides  a  simple, 
fast  and  effective  way  of  generating  low-thrust  trajectories. 

The  procedure  was  found  to  converge  rapidly  for  both  CSI  and  VSI 
vehicles. 

7.  3    The  Mission 

The  mission  selected  was  a  150- day  heliocentric  transfer  from 
Earth  to  Mars.     Flight  time  of  150  days  was  chosen  in  order  that 

numerical  results  could  be  compared  with  a  162  -day  coplanar  trajectory 

12 
generated  by  Friedlander      .     The  period  chosen  corresponds  to  the 

37 
favorable  opposition  of  Mars  during  the  summer  of  1971      .     The  depar- 
ture date,  from  a  position  on  Earth's  sphere  of  influence,  is  J.  D.    244- 
1090.  5.     Arrival  at  Mars  sphere  of  influence  is  150  days  later.     The 
vehicle  is  assumed  to  have  Earth's  orbital  velocity  at  departure  and  to 
match  the  Martian  velocity  at  arrival.     These  velocity  conditions  are 
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not  optimal  for  the  entire  mission,    i.  e.    escape,  transfer,  capture, 
however  the  author  wanted  to  test  out-of-plane  components  of  the  state 
vector.     The  use  of  planetary  velocities  for  initial  and  final  velocity 
permitted  this. 

Both  CSI  and  VSI  modes  of  control  were  tested.     The  assumption 
was  made  that  a  space  vehicle  of  given  size  and  power  could  be  con- 
trolled under  either  mode  of  operation,  CSI  or  VSI.     Thus  two  complete 
sets  of  data  were  generated  which  differ  only  in  mode  of  control.     In 
addition  the  adjoint  sets  for  an  optimal  thrust  program  and  an  optimal 
acceleration  program  were  generated  for  each  control  mode.     The  pur- 
pose of  having  both  types  of  adjoint  functions  is  to  permit  all  engine 
anomalies  of  interest  to  be  studied.     The  subtle  difference  between  the 
two  sets  of  adjoint  functions  is  discussed  in  Appendix  C 

7.  4   Computational  Coordinates 

The  coordinate  system  used  for  the  mission  analysis  is  a  simple 
but  effective  system  which  is  defined  by  the  transfer  plane.     The  trans- 
fer plane  is  the  plane  which  passes  through  the  sun  and  contains  the 
desired  departure  and  arrival  points.     The  x  axis  passes  through  the 
launch  point,  the  z  axis  is  normal  to  the  plane  in  the  northerly  direction, 
and  the  y  axis  completes  the  right  hand  triad.     A  procedure  for  trans- 
forming ephemeris  data  into  computational  coordinates  is  derived  in 
Appendix  F„     Figure  7 -a  illustrates  the  geometry  of  the  transfer. 

7.  5   Engine  Selection 

Engine  sizing  for  the  sample  mission  was  computed  on  the  basis  of 
the  mass  distribution  for  maximum  payload  derived  in  Chapter  II       That 
is 

JV    =    P2£L   SJX   -j  (7-i) 

a  m0  V   a 

A  trial  trajectory  for  a  mass  independent  VSI  rocket  was  computed  to 
obtain  a  first  approximation  for  J       The  value  for  the  one-way  transfer 
was  scaled  up  by  a  factor  of  4  to  provide  a  more  realistic  value  for  a 
round  trip  mission;  also  to  insure  that  mass  would  never  be  reduced  to 
zero  during  computer  tests.     A  value  for  a  of  10  kg/kw  was  selected  as 
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being  reasonable,  though  perhaps  unattainable  in  less  than  a  decade. 
The  resulting  value  for  exhaust  power  is 

_P_    =  o.  0242   —  (7-2) 

m0  kg 

This  value  for  power  was  used  for  both  modes  of  control,  (VSI)  and 
(CSI). 

The  VSI  trajectory  was  generated  and  the  resulting  value  for  initial 
acceleration,   1„  57  X  10       g  ,  was  used  in  the  field-free  space  equa- 
tions of  Appendix  A  to  obtain  an  optimum  initial  acceleration  (thrust 

_4 
level)  for  the  CSI  vehicle.     A  value  of  1.  2  X  10       g    was  used.     This 

value  of  optimum  initial  acceleration,  based  on  field-free  space,  is  in 
surprisingly  good  agreement  with  the  more  rigorous  computations  of 
Melbourne  and  Sauer, 

7.  6    Numerical  Results 


The  computer  output  data  for  the  final  iteration  of  each  mode  of 
control  are  reproduced  as  Appendix  rL     In  addition,  plots  of  several 
interesting  output  quantities  are  presented  as  Figures  H-a  through  H-r. 
The  data  samples  and  the  plots  are  strictly  valid  only  for  the  particular 
cases  they  represent,  however  they  are  indicative  of  the  order  of  mag- 
nitudes applicable  to  many  one-way  missions  in  the  solar  system.     The 
adjoint  functions,  for  example,  (Figures  H-d  through  H-r)  are  in  good 
agreement  with  the  values  obtained  by  Friedlander  for  a  coplanar  trans- 
fer. 

It  is  interesting  to  observe  that  the  acceleration  magnitude  for  a 
VSI  vehicle.  Figure  H-a,  is  approximately  a  linear  function  of  time  and 
almost  symmetric  with  respect  to  the  midpoint,     Comparison  with  the 
field-free  space  acceleration  program,  Figure  2-c,  which  is  linear  and 
symmetric,  confirms  that  to  first  order  analysis  of  VSI  trajectories  in 
field-free  space  may  be  applied  to  the  gravitational  field. 

This  result  only  confirms  the  work  of  other  investigators  and  was 
anticipated.     The  analysis  for  CSI  vehicles  in  field-free  space  is  not 
as  straightforward  and  requires  a  large  amount  of  work.     The  field- 
free  space  derivations  in  Appendix  A  result  in  a  value  for  the  optimum 
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CSI  thrust  level  which  is  in  good  agreement  with  the  computer  results 
of  Melbourne  and  Sauer,     The  value  obtained  for  the  final  mass  ratio 
is  larger  than  the  value  which  results  from  computer  studies  in  the 
gravitational  field.     Referring  to  Figure  2-d,  observe  that  the  point 
representing  the  CSI  transfer  falls  slightly  below  the  curve  represent- 
ing field-free  space  prediction.     This  is  satisfying  in  that  the  penalty 
for  CSI  control  is  less  than  expected.     However,  before  concluding  that 
field-free  space  is  always  a  good  predictor  for  CSI  transfers,  it  will 
be  necessary  to  test  a  large  number  of  additional  points  to  determine 
if  agreement  is  sufficiently  good  to  warrant  the  computational  effort. 
The  limitation  of  time  has  prevented  the  author  from  undertaking  such 
a  study. 

In  Figures  H-b  and  H-c  the  transfer  plane  components  of  position 
and  acceleration  are  plotted.     The  out-of-plane  components  are  quite 
small  and  are  not  considered  further.     A  surprising  result  is  that  the 
maximum  difference  between  the  physical  paths  of  the  VSI  and  CSI 
trajectories  is  less  than  0.  01  A  u.     The  large  difference  in  the  form 
of  the  acceleration  programs,  Figure  H-a,  leads  one  to  anticipate 
rather  large  trajectory  differences.     This  is  found  not  to  be  the  case. 

A  very  important  characteristic  of  continuous -thrust  rockets  is 
shown  in  Figure  H-p.      This  is  the  sensitivity  of  position  error  to  mass 
change  when  the  thrust  program  is  specified.     The  y  component  of 
position  at  arrival  time  is  subject  to  an  error  of  0.  81  X  10     km  for 
1%  variation  in  launch  weight.     For  a  one -ton  vehicle  this  is  approxi- 
mately equivalent  to  a  25,000  mile  terminal  error  for  each  pound  of 
launch  weight  variation.     The  large  sensitivity  to  mass  changes 

emphasizes  the  requirement  for  accelerometers  with  sensitivities  of 

-5 

the  order  of  10       g.     Such  instruments  will  be  extremely  important  in 

the  navigation  of  low-thrust  spacecraft. 

7.  7    Convergence  of  the  Computation  Routine 

The  results  discussed  in  the  preceding  section,  although  interest- 
ing, are  only  by-products  of  the  experiment  devised  to  test  linear 
guidance  theory  as  a  computational,  method  for  trajectories.     Conver- 
gence of  the  routine  is  crucial  to  the  thesis  since  the  linear  theory 
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purports  to  null  the  miss  vector  for  small  variations  in  state.     Appli- 
cation of  the  guidance  theory  to  the  launch  state,  i.  e„    trajectory  com- 
putation, provides  the  most  severe  test  of  guidance.     The  results  of 
this  test  are  plotted  in  Figure  7-b.     The  criterion  chosen  as  indicative 
of  convergence  is  the  RMS  value  of  miss  vector  length.     In  order  to  be 
dimensionally  compatable.   the  velocity  components  are  multiplied  by 
flight  time,   150  days.    Therefore  the  ordinate  in  Figure  7-b  is  the  RMS 

/    ir    \ 

value  in  A.  u.  ,  of  <         >.      >  .       The  abscissa  is  the  iteration  number. 

If  ivJ 
Several  cases  were  tested;  the  three  plotted  are  indicative  of  all  tests. 

The  first  one  or  two  points  in  each  case  represent  the  result  of 
estimating  initial  values  for  the  Euler  variables,,    v_s  and  were  widely 
scattered.     These  points  are  not  plotted. 

The  unrestricted  VS1  trajectory  converges  very  rapidly;  reducing 
|  |  I    by  approximately  two  orders  of  magnitude  each  iteration.     Con- 
vergence is  very  smooth.     The  CSI  cases  converge  less  rapidly  and 
less  smoothly,  but  satisfactorily.     In  general,  as  the  thrust,  limit,  a    , 
is  reduced,  convergence  is  slower.     If,  for  example,  a     is  reduced  e 
below  the  value  which  permits  no  coast  period,  the  procedure  will  not 
converge.     This  physically  corresponds  to  a  vehicle  "in  extremis" 
such  that  insufficient  thrust  is  available  to  complete  the  mission      The 
result  also  indicates  that  reserve  power  must  be  available  for  guidance. 
A  CSI  vehicle  which  is  beyond  the  coast  phase  cannot  correct  for  state 
error  in  all  six  components  of  ^  unless  reserve  power  is  available.    If 

reserve  power  is  not  available,,  guidance  must  be  based  on  a  formula- 

14 
tion  such  as  Pfeiffer's      ,  which  gives  the  minimum  miss. 

It  should  be  reported  that  the  CSI  miss  vector  magnitude  entered 
a  small  limit  cycle  around  the  target  prior  to  addition  of  the  term 
3t,  l$v_.     This  term  is  derived  in  section  6.  3.     Prediction  of  the  switch- 
ing time  by  including  the  above  deriva+ive  is  necessary  to  obtain  the 
desired  accuracy. 

It  may  be  desirable  to  include  an  additional  term  in  the  procedure 
to  smooth  convergence  of  the  CSI  routine.     The  term  9m/  dv,  evaluated 
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Fig.  7-b.    Convergence  of  computation  routine. 
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at  the  switch  points,  will  probably  be  sufficient  to  do  this.  Computer 
tests  were  completed  before  the  significance  of  this  term  was  recog- 
nized and  time  limitation  precluded  retesting  of  the  procedure. 

The  accuracy  obtained  from  the  tests  is  significant.     The  routine 
was  stopped  after  a  predetermined  number  of  iterations  to  conserve 
computer  time,  however  for  VS1  trajectories,  the  largest  error  in  a 
component  of  position  was  of  the  order  of  one  mile.  The  largest  error 
in  velocity  was  of  the  order  of  20  feet  per  hour ■„    These  values  are  given 
in  the  computer  results  in  Appendix  H       The  computer  was  stopped 
before  attaining  this  accuracy  for  CSI  trajectories,  however  the  con- 
vergence plots  in  Figure  7-b  indicate  such  accuracy  is  attainable  with 
a  sufficient  number  of  iterations. 
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CHAPTER  VIII 


SUMMARY  AND  CONCLUSIONS 


8.  1   Summary 

In  this  report  the  objective  has  been  to  justify  the  argument:  "There 
exists  a  linear  method  which  produces  a  prcpellan'1:- optimal  control  pro- 
gram in  a  noniterative  form  for  guidance  of  low-thrust  space  vehicles 
and  which  provides  a  simple,,  rapidly  converging  iterative  technique  for 
computing  propellant-optimal  trajectories,  " 

In  pursuit  of  this  objective  it  is  necessary,  in  Chapter  II,  to  develop 
the  parameters  which  characterize  low-thrust  vehicles.     Due  to  "state 
of  the  art"  restrictions  on  variable -specific -impulse  (VSI)  machines, 
the  parameters  for  constant-specific -impulse  (CSI)  vehicles  are  devel- 
oped in  addition  to  the  idealized  mode  of  contr<  ' 

As  an  adjunct  of  Chapter  II,  the  equations  for  field-free  space  appli- 

-. 
cation  of  CSI  control  are  derived  in  Appendix  A,     This  derivation  is 

used  to  test  the  validity  of  approximating  CSI  transfers  in  the  solar 

system  by  field-free  space  analysis,  as  has  previously  been  done  for 

constant -power  vehicles. 

The  parametric  derivations  provide  a  starting  point  for  investiga- 
ting propellant-optimal  guidance  of  low-thrust  vehicles.     For  this  study 
the  vehicle  is  characterized  by  its  "state".     The  state  is  represented  by 
a  vector  consisting  of  three  components  of  position,  three  components 
of  velocity  and  the  vehicle  mass.     The  differential  equation  of  state  is 
linearized  by  considering  variations  of  the  state  vector  with  respect  to 
an  optimal  reference  trajectory.    Auxiliary  functions,  called,  adjoint 
variables,  are  used  to  solve  the  linearized  diff<  i  <  ntial  equation.     The 
solution  to  the  expanded  equations  is  called  the  state  transition 

equation  or  fundamental  guidance  equation  and  serves  as  the  basis  of 
the  guidance  theory. 
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The  linear  guidance  theory  is  developed  in  Chapter  III  and  applied 
to  both  modes  of  vehicle  control.     The  solution  for  an  unrestricted  VSI 
vehicle  is  derived  using  the  calculus  of  variations  and  then  using  Pontry- 
agin's  maximum  principle  to  prove  that  both  methods  lead  to  the  same 
result.     The  CSI  problem  is  then  solved  using  the  Hamiltonian.     The 
crucial  step  in  either  approach  is  deletion  of  the  adjoint  functions  asso- 
ciated with  the  mass  rate  equation.     Deletion  of  this  set  of  functions  to 
form  a  reduced  adjoint  set  permits  a  solution  to  be  obtained  directly. 
Otherwise,  the  system  is  indeterminate  due  to  a  singular  matrix.     Proof 
of  the  singularity  and  justification  for  deleting  the  one  set  of  adjoint 
functions  are  treated  in  Appendix  C  and  Chapter  V  respectively. 

As  an  adjunct  of  Chapter  III,  explanation  of  the  adjoint  relationship 
and  derivation  of  Pontryagin's  principle  are  presented  in  Appendices  B 
and  D  respectively.     Also  in  Appendix  D  is  a  discussion  of  optimality 
criteria  as  derived  from  the  Pontryagin  principle  and  from  the  calculus 
of  variations.     Important  properties  of  the  state  transition  equation  are 
presented  and  discussed  in  Appendix  C 

The  problem  of  estimating  the  vehicle  state  is  studied  in  Chapter  IV 
as  a  problem  in  navigation.     The  method  of  redundant  measurements 
studied  by  Battin,  Stern  and  Potter  is  extended  to  include  continuous 
measurement  of  low-thrust  engine  performance.     Based  on  the  concept 
of  filtering  redundant  data  with  a  biased  filter,  two  methods  are  derived 
for  predicting  future  state.     One  is  a  simple  method  for  short  term 
prediction;  the  other  is  more  complex  but  also  more  accurate  for  long 
range  prediction.     The  effect  of  omitting  engine  measurements  is  dis- 
cussed. 

The  problem  of  computing  optimal  reference  trajectories  is  dis- 
cussed in  Chapters  V  and  VI.     By  first  deriving  the  Euler  equations  with 
the  calculus  of  variations  and  then  demonstrating  that  the  set  contains 
a  superfluous  Euler  equation,  the  linear  guidance  theory  is  shown  to  be 
useful  for  computing  optimal  reference  trajectories.     The  superfluous 
Euler  variable  is  a  result  of  optimizing  a  state  variable. 
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To  provide  a  comparison  of  linear  guidance  with  better  known  com- 
putation techniques,  the  steepest  ascent  formulation  of  the  trajectory 
problem  is  presented  in  Appendix  E 

Chapter  VII  describes  the  sample  mission  that  was  simulated  on  a 
digital  computer  to  test  the  guidance  equation  as  a  computational  device 
for  trajectories.     The  results  of  the  test  are  discussed. 

8.  2    Conclusions 


On  the  basis  of  the  analysis  and  the  subsequent  simulation    the 
following  conclusions  were  reached. 

1.  The  linear  guidance  method  is  applicable  to  guidance  of  low- 
thrust  vehicles  in  interplanetary  flight. 

2.  The  rapid  convergence  of  the  test  procedure  proves  the  method 
to  be  applicable  for  computing  propellant -optimal  trajectories. 

3        For  CSI  transfers,  reserve  power  for  guidance  is  necessary. 

4.  Mathematical  descriptions  of  many  optimization  problems 
which  contain  a  state  variable  in  the  cost,  will  produce  a 
superfluous  Euler  variable  when  the  classical  calculus  of 
variations  is  used.     Deletion  of  this  Euler  variable  removes 
several  difficulties  associated  with  solving  optimal  control 
problems. 

5.  Sightings  on  celestial  bodies  at  discrete  intervals  may  be 
combined  with  continuous  measurements  of  engine  perform- 
ance to  estimate  the  state  of  the  vehicle. 

6.  The  uncertainty  in  state  is  increased  if  engine  measurements 
are  omitted. 

7.  Field-free  space  analysis  of  CSI.  transfers  provide  a  reason- 
able approximation  for  the  results  in  a  gravitational  field, 
however  additional  study  is  needed  to  ascertain  the  general 
applicability  of  the  method. 

8.  3    Contributions  of  the  Investigation 

The  items  in  the  report  which  are  believed  to  be  novel  are  dis- 
cussed in  this  section. 
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Derivation  of  the  propellant -optimal  control  law  by  the  method  of 
Chapter  III  is  thought  to  be  original.     The  crucial  step  in  the  derivation 
is  the  formation  of  a  reduced  adjoint  set.     This  permits  inversion  of  a 
matrix  which  otherwise  would  be  singular.     Other  approaches  to  this 
problem  usually  require  the  addition  of  artificial  constraints  or  weight- 
ing matrices  which  remove  the  singularity  but  also  change  the  char- 
acter of  the  cost  function.     The  method  of  this  thesis  preserves  the 
cost  function  and  produces  the  optimal  control  directly. 

The  explanation  of  the  singularity  on  the  basis  of  a  superfluous 
Euler  variable  is  thought  to  be  novel.     The  simplification  resulting 
from  deletion  of  the  extra  variable,  by  using  the  reduced  adjoint  set, 
may  resolve  many  difficulties  associated  with  optimal  control  problems 
other  than  the  one  in  this  report. 

The  equations  which  combine  the  celestial  sights  and  the  continuous 
measurements  of  engine  performance  are  not  believed  to  have  been  pre- 
viously derived. 

Finally,  the  method  of  generating  the  optimal  trajectory  is  believed 
to  be  simpler  than  methods  previously  used.     The  method  results  from 
linearizing  the  state  equations  but  using  the  complete  nonlinear  cost 
function. 

8.  4   Recommendations  for  Further  Study 

The  current  research  has  revealed  some  areas  where  additional 
work  might  produce  fruitful  results  and  other  areas  where  information 
is  lacking  with  respect  to  low-thrust  transfers. 

The  computational  method  needs  to  be  expanded  to  the  many  body 
problem.     Preliminary  work  does  not  show  any  monumental  difficulties 
associated  with  such  an  effort.     The  primary  problem  is  to  include  an 
ephemeris  in  the  routine. 

It  appears  possible  to  use  the  method  as  a  search  routine  for  finding 
families  of  trajectories,  with  the  goal  of  defining  optimal  launch  times 
for  low-thrust  missions. 
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A  study  of  methods  to  join  different  segments  of  an  optimal  tra- 
jectory should  be  interesting.     In  particular  the  idea  of  making  only 
one  coordinate  change  between  launch  and  capture  is  appealing. 

Research  to  define  rigorously  the  concept  of  superfluous  Euler 
variables  is  needed.     It  is  doubtful  that  the  singularities  in  all  optimi- 
zation problems  can  be  removed  by  the  method  of  this  thesis,  but  per- 
haps some  can.     Knowledge  of  the  general  applicability  is  needed. 

When  data  on  the  reliability  of  low-thrust  power  plants  and  thrusters 
is  available,  a  study  of  the  effects  of  engine  anomalies  on  the  probability 
of  mission  success  should  be  made. 

A  study  of  VTA  guidance  using  the  "critical"  plane  associated  with 

31 

ballistic  guidance  is  needed      .     This  concept  is  difficult  to  visualize 

in  the  six  dimensional  phase  space  associated  with  the  equations  of 
motion.     However,  the  existence  of  a  concept  analogous  to  the  "critical 
plane"  idea  might  provide  additional  insight  to  the  guidance  problem. 
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APPENDIX  A 
TRANSFERS  IN  FIELD -FREE  SPACE 

A.  1  Introduction 

In  this  appendix  the  equations  for  point-to-point  transfer  of  low- 
thrust  rockets  in  field-free  space  (FFS)  are  presented.     The  optimal 
transfer  problem  for  a  constant  power  VSI  rocket  is  a  straight  forward 
application  of  the  calculus  of  variations.     Examples  of  this  transfer  are 
found  frequently  in  the  literature.     It  is  reproduced  here  for  complete- 
ness. 

Point-to-point  transfer  of  a  CSI  rocket  is  more  difficult  to  compute., 
The  equations  are  derived  in  this  appendix.     Results  of  the  derivation 
were  applied  to  the  sample  mission  in  the  thesis  to  ascertain  the  validity 
of  using  FFS  predictions  in  the  gravitation  field.     FFS  analysis  is  found 
to  furnish  a  reasonable  approximation  for  CSI  transfers,  however  the 
analysis  is  quite  tedious. 

The  author  is  indebted  to  Mr,,    Neal  Carlson  for  checking  the  deriva- 
tion and  for  suggesting  different  methods  of  approach. 

A.  2  Constant  Power  Transfer 

Assume  that  in  FFS  we  desire  to  traverse  the  distance  L  in  the 
time  T  such  that  the  vehicle  begins  and  ends  at  rest  and  such  that  the 
acceleration  integral  is  a  minimum. 

That  is 

T 
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dt  (A-l) 
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Forming  a  functional  F,  one  obtains 

T  T 

'2 

F  -  J   —    dt  +  ?r  (L  -  J     v  dt)  (A-4) 

0        2  0 

Direct  application  of  the  variational  technique  produces  the  optimal 
velocity  schedule  from  which  the  acceleration  is  easily  derived. 

a=%-    (1  -—  )  (A-5) 

T  T 

Evaluating  the  acceleration  integral,  obtain 

2 

;3 


J  =  ^=U-  (A-6 


T 
From  equation  (A-5)  the  optimal  initial  acceleration  is 

ao  ■  P"  <A"7> 

A  plot  of  equation  (A-5)  is  presented  in  the  discussion  of  Chapter  II  as 
Figure  2-c. 

A.  3  Constant- Specific-Impulse  Transfer 

The  well  known  equations  for  a  conventional  rocket  must  be  inte- 
grated in  the  analysis  of  a  CSI  transfer.     That  is 

v  =  c  In  MR  (A-8) 

Since  the  variational  approach  is  degenerate  for  this  type  of  optimiza- 
tion problem,  -we  shall  simply  apply  the  known  result  from  the  Hamil- 
tonian  approach  of  Chapter  III,  that  thrust  is  either  full  on  or  full  off  for 
CSI  vehicles  in  a  linear  field.     Integration  of  equation  (A-8)  is  over  the 
initial  and  final  thrusting  periods  with  a  coast  in  between.     That  is 

ti  T 

L  =  /  v  (t)  dt  +  /    v  (t)  dt  +  v  (t   )  (t     ■-  t   )  (A-9) 

0  t2 
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subject  to  the  boundary  conditions  (A-3).      This  rather  lengthy  compu- 
tation reduces  to 

L  =  c  (t2  -  tx)  In    I— 4 j     +aQt12  (A-10) 

where  a     in  the  initial  acceleration  and  m  is  the  normalized  flow   rate 
o 

(i.e.    m     =  1).     Some  additional  algebraic  manipulation  allows  (A-10) 
to  be  cast  in  the  form 

-2^-  =  -   [mT  +  m2  ~  l\  In  rr^  +  (1  -  m^2  (A-ll) 


where 


m1  =  1  -  mt1  (A-12) 


From  the  work  of  Chapter  II  observe  that  for  m  a  positive  number 

a     =  mc  (A-13 

o 

p  =•+-  mc2  (A- 14) 

2 

Thus  2 

a 
m  =—2—  (A-15) 

2p 

c  =  ^  (A-16) 

a 
o 

Equation  (A-Tl)  may  now  be  written  in  terms  of  initial  acceleration, 
power,  mass,  and  the  dimensions  of  the  problem,  L  and  T„ 


a3  L  /a2  T2 

o  /o  .       2       A      -,  ,   ,  1  ,2 

— 2~    =  -         — —    +  m     -  lj      In  nij  +  (1  -m^) 


... 

4p  V     2p  /  (A-17) 


fi  T 

In  the  VSI  case  a     was  — s-  .     Let  the  CSI  initial  acceleration  be  an 

O  rp2 

unknown  multiple  of  that  value 

aQ  =  x^  (A-18) 
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Further,  define  a  parameter  R. 

2 

.3 


1  RT  ' 

R  i    -^V  (A-19) 


^-R     x     =  -(  R  x     +  m     -  lj  In  m     +  (1  -  m..)' 


pT 

Substituting  (A- 18)  and  (A-19)  into  (A- 17)  one  obtains 

(A-20) 

From  the  condition  that  the  velocity  change  must  be  the  same  for  both 
thrusting  intervals  if  boundary  conditions  are  satisfied,  the  final  mass 
ratio  is 

1       2 
MR=(^-)  (A-21) 

Therefore  if  m     is  a  maximum,  mass  ratio  is  a  minimum  and  likewise 

for  the  acceleration  integral. 

3m 

Using  equation  (A-20),  solve  for  x  such  that —    =  0.     The  result  is 

dx 

3m. 

=  0  (A-22) 

dx 

when 

Rx(JRx+2  1n  mx)  =  0  (A-23) 

x  =  0  is  a  trivial  solution  and  the  desired  solution  is 


x  =  -  -    In  m,  (A-24) 

R  l 


x  =   -  In  MR  (A-25) 

R 


If  equation  (A-24)  is  substituted  into  (A-20)  the  result  is 

-  —  (In  m   )3  +  (1  -  m2)    In  m     +  (1  -  m2)  =  0 
3R 


(A-26) 


The  transcendental  equation  (A-26)  is  the  relation  between  m..   and  R 
for  optimal  CSI  operations  in  field-free  space.     The  optimum  initial 
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acceleration  is  given  by  (A-24)  or  (A-25).     Equation  (A-26)  is  plotted 
in  Figure  2-cL     It  is  observed  that  CSI  control  always  uses  more  pro- 
pellant  than  VSI  control.     Also  plotted  in  Figure  2-d  is  the  CSI  curve 

fi  T 

for  a     =  — r~  (i.  e.    x  =  1),     This  value  of  a     is  considerably  more  expen- 

O         T  O 

sive  than  a  transfer  with  a     =  a     .      ,.  .     Figure  A-a  contains  plots  of 

o         o  (opt)  ^ 

the  required  coast  time  for  CSI  control.     These  are  obtained  by  solving 
equation  (A- 10)  for  (t„  -  t-.  )/T  in  terms  of  m..   and  R.     Again  the  work 
involved  is  substantial. 

Observe  that  for  the  case  where 

m1  =  1  -e  (A-27) 

where  e    is  a  small  quantity,,  equation  (A-27)  may  be  substituted  into 
(A-26)  and  the  resulting  expression  solved  for  m..   by  neglecting  higher 
order  terms. 

m     =  i  -  3R_  (A-28) 

1  16 

By  approximating  lnm.   as 

In  mi  =  -  e    (1  +  1/2  e    +••••••)  (A-29) 

one  obtains 

x      ,  =  3/4  (A-30) 

opt         ' 

Therefore,  for  transfers  such  that  the  propellant  consumption  is 
small  with  respect  to  the  total  initial  mass,  the  optimum  initial  acceler- 
ation for  CSI  vehicles  is  approximately  three -fourths  of  the  correspond- 
ing optimum  initial  acceleration  for  VSI  vehicles. 
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APPENDIX  B 
ADJOINT  SYSTEM  OF  DIFFERENTIAL  EQUATIONS 

B„  1  Introduction 

The  purpose  of  this  appendix  is  to  review  briefly  the  concept  of  the 
adjoint  system  of  differential  equations  and  to  illustrate  its  use  in  a 
simple  guidance  problem.     Although  the  concept  is  not  new,  it  has  not 
been  used  as  a  standard  engineering  technique  in  guidance  and  control 
problems  until  very  recently.     For  solution  of  the  two-point  boundary 
value  problem  which  occurs  in  the  navigation  and  guidance  of  space- 
craft however,  the  method  of  adjoints  is  a  particularly  useful  tool. 

The  development  in  this  section  is  not  intended  to  be  rigorous  and 
exhaustive  but  merely  illustrative  of  the  method  used  throughout  the 
thesis.     The  material  is  taken  primarily  from  the  lecture  notes  of 
Professor  Frank  D.    Faulkner  of  the  United  States  Naval  Postgraduate 
School. 

B.  2  Method  of  Adjoints 

Consider  the  following  ordinary  scalar  differential  equation, 

x  +  3x  +  2x  =  f(t)  (B-l) 

which  may  be  used  to  define  the  operator 

2 

L(x)    =  (-^r+  —    +  2)  x  =  f(t)  (B-2) 

dt  dt 

where  f(t)  is  an  arbitrary  but  known  function.     In  a  manner  closely  re- 
lated to  the  method  of  variation  of  parameters,  form  the  integral 

T  T 

J  =    f     ALIx)  dt    -     J     Af(t)  dt  (B-3) 

0  0 

where  TV   is  an  unspecified  function  which  will  be  chosen  later  to  satisfy 
certain  boundary  conditions  in  addition  to  a  functional  relationship.     In- 
tegrating equation  (B-3)  by  parts  to  eliminate  the  derivative  of  x  from 
the  integrand,  one  obtains 


101 


T         T 

J  =  (  Ax+  3  Ax  -    Ax)   |      +/    x(A-3A+2A)dt 

0  (B-4) 

If  J\_  is  chosen  such  that  it  satisfies 

L*  (  A)  =  (-^-    -    —    +  2)    A=  0  (B-5) 

dt  dt 

then  the  definite  integral  J  is  a  function  only  of  the  first  term  on  the 
right  side  of  equation  (B-4).     The  operation  of  integrating  by  parts  to 
eliminate  the  dependent  variable  from  the  integrand  defines  the  adjoint 
operator  L*.     If  L*  -  L  the  system  is  said  to  be  self  adjoint  and  has 
some  very  useful  properties.     However,  these  are  of  no  concern  at  the 
moment. 

A  general  solution  to  the  adjoint  equation  (B-5)  is 

A=  Cx  e1  +  C2  e2t  .  (B-6) 

Suppose  that  the  desired  quantities  are  x(T)  and  x(T)  and  that  x(0)  and 
x(0)  (i.  e.    two  constants  of  integration  associated  with  the  original  dif- 
ferential equation)  are  known.     Equation  (B-4)  may  be  rewritten  as 

T 
[Ax  +-  (3 A  -  A)x]  =  [  Ax+  (3  A-  A)x]  +    /  Af(t)  dt 

t=T  t  =  0     0  (B-7) 

If  one  specifies  that 

A(T)    =  1  (B-8) 

3  A(T)     -     A(T)    -  0  (B-9) 

then  equation  (B-7)  gives  for  x(T) 

T 

x(T)    =  [  Ax  +  (3A  -  A)  x]  +    /   Af(t)  dt  (B-10) 

t=0     0 

The  constants  in  equation  (B-6)  may  be  evaluated  from  the  boundary 
conditions  of  the  adjoint  variable,  equations  (B-8)  and  (B-9),     Designate 
this  solution  A  ■.  .      By  performing  the  algebra  one  obtains 

A1  =  2e2<t"T)    -e'-T  (B-ll) 
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Equation  (B-10)  may  now  be  solved  explicitly  for  x(T). 

T 
x(T)  =  (2e"2T  -  e~T)  x     +  2(e~2T  -  e~T)  x    +  /  A«f(t)  dt 

U  (B-12) 

To  obtain  x(T)  it  is  necessary  to  specify  different  boundary  condi- 
tions for   A.     Designate  the  solution  satisfying  these  conditions  as  A  „. 

A2(T)    =  0  (B-13) 

3  A2(T)    -    A2(T)     =   1  (B-14) 

By  solving  equation  (B-6)  subject  to  the  boundary  conditions  of  (B-13) 
and  (B-14)  one  obtains 

A      ^et-T-e^t-T) 

2  (B-15) 

Inserting  this  result  into  equation  (B-7)  yields  for  x(T) 

T 
x(T)  =  (e~T   -  e~2T)  xQ  +  (2e~T  -  e"2T)  xQ  +  /A2f(t)  dt 

0  (B-16) 

The  preceding  development  illustrates  the  general  technique  to  be 
used.     However,  in  order  to  more  closely  approach  the  formulations 
used  in  the  thesis,,  it  is  convenient  to  write  the  basic  second  order  dif- 
ferential equation,  (B-l),  as  two  first  order  equations. 

x  =  y  (B-17) 

y  =  -  3y  -  2x  +  f(t)  (B-18) 

Associate  with  the  two  equations  two  sets  of  adjoint  functions  A  . ,  and 

A  .„  where  the  subscripts  1  and  2  refer  to  the  associated  equation  and 

the  subscript  i  will  refer  to  the  boundary  conditions  which  will  eventually 

be  assigned  to  the    A's.      The  integral  J  may  now  be  formed. 

T  T 

J  =  /  An  (x  -  y)  +  A      (y  +  3y  +  2x)  dt  =   /  A      f(t)  dt 

0  l  0,  1Z  (B-19) 

Integrate  by  parts  to  obtain 

T        T 

J=(Anx  +  Ai2y)   l0-o/  x(Ail-2yVi2)+y(Ai2-3Ai2  +  Ail)dt 

(B-20) 
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In  order  to  eliminate  the  integral  term  it  is  clear  that  the  adjoint  vari 
ables  must  satisfy  the  relationship 


A 


il 


2  A,  9  ;  ° 


i2 


Ai2-3Al2+An  =  o 


(B-21) 
(B-22) 


By  reducing  the  two  equations  to  one  second  order  equation  and  solving 
subject  to  the  boundary  conditions 


and 


An  (T)  =   1 

A12  (T)  =  0 

A21  (T)  -  0 

A22  (T)  =   1 


(B-23) 
fB-24) 

(B-25) 
(B-26) 


it  is  apparent  that  the  adjoint  system  of  equations  (B-21)  and  (B-22)  is 
identical  with  that  obtained  by  the  first  formulation  of  the  problem. 

In  general,  if  any  set  of  first  order  differential  equations  may  be 
written  in  the  form: 


Xl 


x„ 


>    - 


X 

\     n  J 


H 


x 
v.    n 


fx(t) 


f   (t) 
^  n    V 


(B-27) 


or  equivalently 


x  -  Ax  =  f(t) 


(B-28) 


Then  there  exists  a  set  of  adjoint  functions  satisfying  the  relations 


A  4  AA  =  O 


n 


A(T)    =  I 


such  that 


n 


T 


x(T)    =  A(tx)  xftj)    +  /    A(t)f(t)  dt 


(B-29) 
(B-30) 

(B-31) 
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where  _A_  (t)  is  the  n  by  n  matrix  of  adjoint  functions,  A  is  the  n  by  n 

matrix  of  coefficients  in  the  original  set  of  differential  equations  and 

I     is  the  n  by  n  identity  matrix, 
n  J 

Unfortunately,  in  most  problems  of  interest  the  adjoint  functions 

cannot    be  obtained  in  closed  form  and  it  is  necessary  to  integrate  the 

2 
n     equations  of  (B-29)  backwards  in  time  from  T  to  t.  numerically, 

using  (B-30)  as  the  initial  condition. 

B.  3  Use  of  the  Adjoint  Method  in  a  Simple  Guidance  Problem 

In  this  section  the  set  of  adjoint  functions  for  a  power-limited 
space  vehicle  in  two  dimensional  field- free  space  will  be  developed 
analytically.     The  purpose  is  twofold:     1)  to  further  illustrate  the  tech- 
nique of  section  B.  2,  and  2)  to  provide  a  simple  analytic  example  for 
use  in  illustrating  the  guidance  methods  expounded  in  the  body  of  the 
thesis. 

The  fundamental  guidance  equation  for  thrust- limited  vehicles  has 
been  defined  in  earlier  sections  as 

H 

6sf  =   A,   6s.  +    /   AB  6f  dt  (B-32) 

i  t        t      t 

For  an  actual  interplanetary  transfer  the  A  matrix  will  be  obtained 
simultaneously  with  the  desired  trajectory.     In  any  event,  if  the  tra- 
jectory is  known  the  adjoint  functions  are  also  known  or  may  be  easily 
determined.     To  illustrate  this  assume  a  transfer  in  field-free  space 
such  that  motion  is  along  the  positive  x  axis  commencing  from  rest  at 
the  origin  at  t  =  0.     At  t  =  tR  the  thrust  is  reversed  to  effect  rendezvous 
with  point  xf  at  time  t  =  tf.     See  Figure  B-a. 

y 


t  =  o  t  =  tR  t  =  tf 


^J>  ,xr         CD**    |xf 

L J ^  x 

Figure  B-a 
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The  engine  produces  constant  thrust  along  the  nominal  trajectory  (which 
is  not  optimum).     Presumably  the  engine  must  be  controlled  in  some 
way  to  provide  steering  but  it  is  not  necessary  to  specify  the  method 
of  control  in  order  to  determine  the  TV  matrix. 


The  equations  of  motion  and  constraint  are: 


where 


X   =    V 
X 

y    =    V 

.        y 

v     =  a  =  2p/cm 

v     =  0 

y 

m  =  -2p/c 

p  = 

_  '     2 

=  (exhaust  power) 

(B-33) 
(B-34) 
(B-35) 
(B-36) 
(B-37) 

(B-38) 


2 

c  =  exhaust  velocity 

Since  the  velocity  and  position  along  the  nominal  trajectory  are  easily 
computed  as  functions  of  time  by  direct  integration,  they  may  be  con- 
sidered as  known  functions  and  we  may  proceed  directly  to  the  varia- 
tional equations  which  are  of  primary  interest.     Taking  the  variations 
of  equations  (B-33)  through  (B-38)  yields 

6x  =  6V  (B-39) 

.X. 

6y  =   6Vy  (B-40) 

6vx  =  -     -^-6m    +    i-    6(~lE-)  (B-41) 

cm  m  c 

5v     =  0  (B-42) 

y 

6m  =  0  -     6(^B  )  (B-43) 

c 

These  may  be  formed  into  the  matrix  equation  (B-28)     , 

6s  =  A  6s  +  6f  (t)  (B-44) 
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or 


6s 


0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0      0 
0 


2p/cm' 


0      0      0      0      0 
L0      0      0      0      0 


(   o 


6s  +    < 


0 

_L 
m 

0 


6(2p/c; 


(B-45) 


I     ~  6(2p/c    ) 


The  column  vector  containing  the  power  and  exhaust  velocity  varia- 
tions corresponds  to  the  arbitrary  forcing  functions,  _f(t),  of  section  B.  2 
and  to  the  product  B  6f  in  the  general  guidance  problem.     They  may  be 
written  in  any  manner  convenient  for  the  investigation  at  hand.     They 
are  of  no  further  concern  at  the  moment  and  will  be  carried  along  as 
B6f_. 

The  adjoint  functions  may  now  be  determined  by  direct  integration 
of  equation  (B-46). 

A+  A  a  =  or 


A(T)  =  i5 

Expanding  the  above  and  integrating  the  twenty- five  equations  one  ob 
tains 


(B-46) 
(B-47) 


A  = 


1         0     (tf-t)  0 


(■ 


m 


m 


f 
m 

2m, 


m     , 

+  In  — —     t  >   tD 
m    /  R 


m, 


m 


m 


+ 


R 


m 


In 


mm 


£-V<t 


(tf-t) 


1_ 
m 

m 


1 
m 


f  / 


1 


R 


m, 


t  >  t 


m 


R 

t  <   t 


R 


(B-48) 

Thus  the  adjoint  set  has  been  determined  for  a  simple  case  and  may 
be  used  in  the  study  of  guidance  techniques. 
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APPENDIX  C 

PROPERTIES  OF  THE  FUNDAMENTAL 
GUIDANCE  EQUATION 

C.  1  Summary 

The  contents  of  this  appendix  include  discussions  relative  to  inter- 
changing the  control  vector  and  relative  to  properties  of  the  adjoint  set. 
In  the  absence  of  control  perturbations  the  adjoint  set  alone  describes 
the  effect  of  state  perturbations  and  thus  may  be  called  the  "state  tran- 
sistion  matrix",     Several  of  its  interesting  properties  are  discussed. 

Proof  of  singularities  in  the  guidance  equation  is  presented  in 
section  C.  5. 

C.  2  General  Remarks 


To  provide  a  better  understanding  of  the  methods  used  in  the  solu- 
tion of  the  equations  of  motion,  and  to  show  the  relationship  between 
trajectory  determination  and  guidance  it  is  beneficial  to  examine  vari- 
ous formulations  of  the  variational  equations  of  motion  and  the  con- 
straining equations, 

A  by-product  of  using  the  adjoint  method  with  the  calculus  of  varia- 
tions to  find  an  optimum  trajectory  is  the  solution  of  the  adjoint  set  of 
equations.     Physically  the  adjoint  set  represents  "sensitivity"  coeffi- 
cients or  "influence"  functions  which  permit  the  investigator  searching 
for  a  trajectory  to  determine  the  changes  in  trajectory  parameters 
which  will  move  the  solution  in  the  direction  of  the  desired  optimum. 
It  is  interesting  to  note  that  along  the  optimum  trajectory  the  Lagrange 
multipliers  are  linear  combinations  of  the  adjoint  variables.     From 
the  viewpoint  of  one  searching  for  a  trajectory,  the  adjoint  set  has 
served  its  purpose  once  the  optimum  trajectory  is  determined,     How- 
ever, for  the  guidance  analyst  the  adjoint  functions  serve  a  most  useful 
purpose  by  showing  the  effect  of  spacecraft  perturbations  on  the  final 
state  of  the  vehicle. 
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C.  3  Interchanging  the  Control  Vector 

The  equations  of  motion  and  the  mass  rate  equations,  expressed 
in  a  nonrotating  frame  with  origin  at  the  central  body,  are 

r  =  v  (C-l) 


v  = 


M 


r  +  a 


(C-2) 


m  =  gm  (m,  a) 


(C-3) 


or 


m  -  gm  (f_) 


(C-4) 


The  variables  m  and  a  have  a  functional  dependence  which  must  be  ex- 
plicitly taken  into  account.     In  the  subsequent  discussion  only  the  VSI 
case  will  be  derived  since  CSI  control  follows  an  analogous  argument. 

The  relationship  between  m  and  a.  may  be  written  in  any  manner 
which  is  convenient  to  the  argument  at  hand.     If  one  is  interested  in  the 
effect  of  acceleration  changes  irrespective  of  their  cause  then  write 


m 


2       2 
a     m 

2p 


(C-5) 


The  analysis  of  Chapter  II  shows  that  equation  (C-5)  is  the  correct 
form  for  mass  rate.     By  taking  the  first  variations  of  equation  (C-l) 
through  (C-3)  the  state  variational  equation  is  obtained 


where 


A  8s  +     B  6a 


a 


"O. 


A 


3 
O 


T 
O1     O 


O 


O 

2 

•a^m 

P 


a 


B  - 


°3 


-m' 
l_~P~ 


(C-6) 


(C-7) 


(C-8) 
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where  the  pre -subscript  denotes  that  a  is  the  controls     The  submatrices 
of  A  and  B  are  the  three  by  three  null  matrix  0„;  the  three  by  three 
identity  I„;  the  three  by  three  symmetric  matrix  of  gravitational  grad- 
ients G;  and  the  three  by  one  null  vector  O.     If  the  adjoint  set  for  equa- 
tion (C-6)  is  formed  such  that 

(C-9) 

(C-10) 


A=  -A  A 

A  (tf )  =  I 

then  the  solution, W .3 

must  have  the  form 

AnA-12 

O 

aA  = 

A21A22 

o 

OT      QT 

A  0 

33 


(C-ll) 


where  the  matrix  is  partitioned  into  three  by  three  submatrices^  three 
by  one  vectors  and  a  scalar,  A 

The  form  becomes  apparent  when  equation  (C-9)  is  expanded  and  an 
initial  condition  is  applied,     Notice  also  that  for  the  form  (C-ll)  to 
exist  it  is  not  critical  thatA(tf)    =  Is  only  that  at  some  time  t,A(t)  =1 

The  expansion  of  (6-9)  yields 


7Vn=A12G 


A 


12 


-A 


11 


A13  =  ° 


A21  =A22  G 


A 


22 


-A 


21 


A23 

=  _0 

A  31 

=As2G 

aJ2 

T 

■-A31 

0 

A   „„ 

a^m    » 

33 


P 


33 


(C-12) 

(C-13) 
(C-14) 
(C-15) 
(C-16) 
(C-17) 
(C-18) 
(C-19) 
(C-20) 
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The  form  of  equation  (C-Tl)  is  assured  if  at  some  time  A  =  I. 


Now  rewrite  the  functional  dependence   of  a  and  m  as 

JL 

a  -  — 
—      m 


(C-21) 


m  = 


i_ 

2p 


where  again  1,  m  and  p  are  per  unit  initial  mass, 
must  be  modified  to 


(C-22) 
The  A  and  B  matrix 


-A  = 


°3 

h 

G 

°3 

oT 

oT 

o 

a 
m 
0 


(C-23) 


f8 


O, 

1 
m     3 


L 


p 


(C-24) 


where  the  subscript  denotes  that  f_  is  the  control,     The  state  variational 
differential  equation  is 

6s  =  fA6s  +  fB6f_  (C-25) 

Clearly  the  optimal  trajectory  must  be  independent  of  the  particular 
control  variation  used  in  the  computation  scheme,     That  is,  if  the  varia- 
tional problem  for  computing  trajectories  is  solved  using  f  as  the  con- 
trol,, it  must  result  in  the  same  trajectory  that  would  be  obtained  if  a 
is  used  as  the  control.,     Thus,  two  different  formulations  are  available 
which  may  be  used  at  the  discretion  of  the  investigator.     Since  the  cost 
has  been  specified  as  a   /2  it  is  more  convenient  to  use  a  as  the  control 
for  computing  trajectories.     Nevertheless,  the  use  of  _f  as  a  control  has 
important  significance. 

The  adjoint  set  which  corresponds  to  using  £  as  the  control  will  now 
be  determined.    Expanding  (C-9)  using(C-23)  one  obtains  differences 
from  the  set  (C-12)  through  (C-20)  only  in  the  last  column  of  the  solu- 
tion VV. .  That  is 
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A  13     A12  m 


A-  23     A22  m 


A33  =  0 


(C-26) 
(C-27) 
(C-28) 


For  boundary  conditions  TV  (tf)  =  I,  the  solution  has  the  form 


•A 


"A  1 1  A  i  o  T\ 


11 


12  -a  13 


A21  A22-A-23 


oT    oT 


1 


(C-29) 


The  solutions     A.  and  f  A  differ  only  in  the  last  column.     The  remaining 
numbers  of  the  array  are  identical  for  each  point  in  time. 

The  physical  significance  of  this  difference  is  apparent  from  the 
following  physical  reasoning.     Consider  the  adjoint  set  as  an  array  of 
partial  derivatives:     If  A.  (tf )  =  I  then  the  last  column  is 


A23) 


A 


33 


(C-30) 


If  a  is  the  control  then  one  should  expect  A  1  o  and  A  00  to  be  null  vec- 
tors, since  for  a  given  acceleration  program  variations  in  mass  can- 
not affect  the  end  point  in  position  and  velocity.     However  the  end  point 
for  mass  must  change  because  it  requires  a  different  amount  of  pro- 
pellant  for  different  size  vehicles. 

If  _f  is  the  control  and  the  thrust  program  is  given,  then  changes  in 
mass  will  affect  the  end  points  in  position  and  velocity  due  to  the  change 
in  acceleration.     For  a  given  thrust  program  the  final  mass,  however, 
can  change  only  if  the  initial  mass  changes,  therefore  A  no  -  1. 
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It  is  concluded  that  the  control  vectors  6f  and  6a  can  be  interchanged 
in  the  guidance  equation  simply  by  changing  the  last  column  of  the  ad- 
joint set  and  changing  the  matrix  B.     Along  the  optimal  trajectory  all 
of  these  functions  are  known  thus  the  interchange  can  be  easily  made  at 
the  discretion  of  the  investigator. 

The  purpose  of  interchanging  control  vectors  is  simply  stated.    Com- 
putation of  trajectories  is  most  conveniently  carried  out  with  a  as  the 
control,  due  to  the  formulation  of  the  cost.     Assume  however,  that  it 
is  desired  to  investigate  the  effects  of  changes  in  engine  performance 
or  to  investigate  variables  other  than  a  as  quantities  to  be  measured. 
Then  it  is  convenient  to  write,  from  Chapter  II 

_f  -  mc  (C-31) 

then  6f  =  m  6c  +  c  6m  (C-32) 


or 

Further,  let 

then 
or 


6f 


6f 


6f_ 

=  [c  mL] 

f  -  "2p 

2 
c 

2c 

2p 

2  c  c 


T 


_2_ 
2 


(C-33) 
(C-34) 

(C-35) 
(C-36) 


If  in  equation  (C-25)  or  in  the  guidance  equation,  6f^  is  replaced  by 
either  (C-33)    or  (C-36)  it  is  possible  to  study  the  effects  of  the  par- 
ticular variations  represented  without  computing  a  new  adjoint  set. 

The  ability  to  interchange  control  vectors  at  will,  without  requir- 
ing a  complete  new  set  of  computations  makes  the  fundamental  guidance 
equation  a  very  powerful  tool. 

C.  4  Properties  of  the  State  Transition  Matrix 

If  the  fundamental  guidance  equation  is  written  for  two  different 

times  t  -  t.  and  t  =  t.  and  subtracted,  then 
i  J 
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or 


A(t.)  6s  .  =  A(t.)  6s  .  +     f  AB  6a  dt 
i-i  J      -J        f 


t. 


6s  .  =  T. .  6s  .  +      A     (t. )     f    A  B  6a  dt 

—  i         n—1  J  v       i     J,  — 


i] 


(C-37) 


(C-38) 


where 


T..  =  A     (t.)  A(t.) 
ij  i  J 


(C-39) 


In  the  absence  of  control  perturbations,  T. .  completely  describes  the 

transformation  of  the  state  vector  at  t.  into  the  state  vector  at  t..     Thus 

J  i 

T. .  is  called  the  "state  transition  matrix".     It  is  necessary  to  show 

ij  J 

that  T.  .  is  nonsingular.     It  is  desirable  to  show  that  it  can  be  inverted  by 

ij  -1 

inspection  and  that  T. .     =  T...     This  last  property  is  readily  observed  by 

J  J 

taking  the  inverse  of  equation  (C-39). 


•V1  =  '^V 


-1 


(C-40) 


(T..)"1  =  T..   =  A  "VV 

ij  Ji  J         i 


(C-41) 


To  examine  the  other  properties,  it  is  necessary  to  look  first  at  A . 
From  section  C.  3  write  A   in  its  general  form 


A: 


A 


11 


A 


12 


A 


13 


■A- 21     ^22      ^23 


O 


T 


Ot    A 


33 


(C-42) 


Equation  (C-42)  is  to  be  understood  as  representing  both  the  adjoint  set 
with  f  as  control  and  the  set  for  a  as  control.  The  last  column  takes  on 
the  value  appropriate  to  the  control. 

Now,  partition  A  and  consider  only  that  portion  defined  by 

An     A12 


H. 

J 


LA21 


A 


22  J 


(C-43) 


.1 
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where  the  subscript  refers  to  time  t..     The  differential  equations  for 
this  subset  of  A  are  reproduced  from  section  C.  3. 


An=A12G 


(C-12) 


A 


12 


A 


A  21   =  A  22 


11 
G 


(C-13) 
(C-15) 


A 


22 


-A 


21 


(C-16) 


Observe  that  these  equations  are  independent  of  the  control  and  are 
completely  specified  by  the  physical  path  and  the.  boundary  conditions. 

Using  C-T3  and  C-165  H.  may  be  rewritten  as 


H.  = 
J 


-A12     A12 


~A22    A 


22 


(C-44) 


Now  define  a  new  matrix  H.: 

J 


H* 

J 


Aj2     A  22 


A 


12 


A 


22  _ 


(C-45) 


H.*  and  H.  are  related  by 
J  J 


H 


,T 


=  H.P 

J 


(C-46) 


where 


"O, 


h 

ex 


(C-47) 


P  is  a  skew  symmetric  matrix  and  possess  the  properties: 


P     =  -I 
PT  -  P"1  = ■  -  P 


(C-48) 

(C-49) 


Since  in  equation  (C-46)  t.  is  an  arbitrary  time  along  the  trajectory, 

J  T 

the  subscript  may  be  dropped  and  H*      considered  as  a  time  varying 

matrix. 
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H*T  =  HP 


(C-50) 


Differentiating: 


"  ,  T 
H*      =  HP 


(C-51) 


From  equations  (C-12)  through  (C-16)  and  (C-44)  it  is  possible  to 
write 


H  +  H  A*  =  O, 


where 


A*  = 


°3         h 


■G        O 


3  J 


(C-52) 


(C-53) 


Substituting  (C-52)  into  (C-51)  and  using  (C-50)  one  obtains 

•  rp  rp 

H*      =  +  H*      P  A*  P 

Multiplying  out  the  known  product  yields 


PA*  P  =  A 


.„T 


(C-54) 


(C-55) 


since  G  is  a  symmetric  matrix.     Thus  the  transpose  of  equation  (C-54) 

is 


H*  =  A*  H* 


(C-56) 


Comparing  (C-52)  and  (C-56)  it  is  observed  that  H  and  H*  satisfy  the 
adjoint -fundamental  relationship  mentioned  in  Chapter  III.     That  is 


dt 


(H  H*)  =  Q6 


or  integrating 


where 


H  *  =  H.  H* 

f  J       J 


H  (tf)  =  Hf  =  Ig 

Substituting  equation  (C-50)  into  (C-58)and  using  (C-49)  obtain 

T  T 

-PH.    =  -  H.  PH 
f  J         J 


(C-57) 


(C-58) 


(C-59) 


(C-60) 
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Using  (C-59) 


H.  PH 

J         J 


T 


(C-61) 


Any  matrix  H  exhibiting  the  property  of  equation  (C-61)  is  said  to 
be  symplectic  and  its  inverse  is  given  by 

1 


H 


T 
PH      P 


(C-62) 


The  inverse  of  the  submatrix  H  of  A  has  been  obtained  by  rearrange 
ment  of  the  components.     To  obtain    A       two  cases  must  be  considered: 
the    adjoint  sets  for  the  control  a  and  for  the  control  f . 


If  A  =      A  ,  then 

ex 


.A 


H 

T 


°6 


°6 

A33 


where  Ofi  is  a  6  component  null  vector.      Further 

-1 


A 


-l 


H 


O; 


°6 


A 


33 


If  A  =  f  A,  then 


A  = 


H 


A 


O 


T 


(C-63) 


(C-64) 


(C-65) 


where   A.  fi  denotes 


Equation  (C-65)  may  be  written  as  the  product  of  two  matrices, 


A  = 


Then 


A" 


_-:6 

A6" 

""-H 

9S~ 

_£ 

1 

L$ 

1 

"-H"1 

°6 

-I 

A6 

T 
_-6 

1 

T 
_^6 

1 

(C-66) 


(C-67) 
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A 


-l 


H"1        -H_1A  6 


o 


(C-68) 


since  the  second  term  in  (C-67)  is  its  own  inverse.     Therefore _A_  has 


now  been  inverted  for  all  cases.     Now  consider  T. 


ij 


t. .  =a:1  a  . 
IJ       1        J 


(C-39) 


For  A    =     A  ,  use  (C-63)  and  (C-64) 

3. 


T.  . 
a    ij 


PL1  H.i       Oc 


1  A 

I  7V33  j 


°6     ,7U 


33  i 


For  A    =  fA,  use  (C-65)  and  (C-68) 


(C-69) 


,T.. 


H.  X  H. 

i         l 


T 
°6 


H._1<;Afi.  -ar 

i    V-6i       -  6j 


(C-70) 


Finally  it  is  necessary  to  show  that  T.  .  is  nonsingular.     To  do  this 


ij 


-1 


it  is  sufficient  to  show  that  the  determinant  of  H.     H.  is  never  zero  and 

i        J 
that  A  oo  is  never  zero. 

From  equation  (C-20)  it  is  apparent  that  if  Aoo  (0)  =  1  it  increases 
monotonically,  thus  is  never  zero. 


Now  examine  the  determinant  of 


1 


H.  ^H. 

i        J 


From  (C-61) 


det  HPH      =  det  P  =  1 


(C-71) 


Since  the  determinant  of  a  product  is  the  product  of  the  determinants 


(det  H)     =   1 


(C-72) 
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From  equation  (C-62) 

det  H"1  =  det  (-PHTP) 
det  H_1  -  det  H 
since  det  (-P)  -  (-l)k  det  P  =  1 
where  k  =  6  is  the  order  of  P. 


(C-73) 
(C-74) 
(C-75) 


Since  the  subscripts  i  and  j  are  completely  arbitrary  and  since  H 
is  continuous, 


det  HT1  H.  =  det  H"1  H  =  1 

i         J 


(C-76) 


Thus  T,  .  is  nonsingular. 


C.  5  Singularities  of  the  Guidance  Equation 

Solution  of  the  guidance  equation  for  a  control  program  which  will 
satisfy  certain  terminal  conditions  invariably  results  in  the  inversion 
of  a  matrix  integral  expression.     In  section  3.  8  the  solution  for  FTA 
guidance  requires  such  an  inversion,  namely 


M 


-1 


tf 


/    A*  BB1  A*     dt 


(C-77) 


In  this  section  the  proof  will  be  given  that  M  is  not  singular  for 
(t~  -  t,)  >  0    but  that  the  corresponding  matrix  denoted  by  M  ,  with  A. 
in  place  of  A*,  is  singular  in   the  case  of  unconstrained  control  pro- 
grams, for  all  values  of  t1  along  the  optimal  trajectory.  We  shall  pro- 
ceed to  prove  that  M     is  singular  and  by  induction  show  that  M  is  not 
singular. 

The  following  proof  of  the  singularity  in  M     was  suggested  by  Dr. 
James  Potter. 

If  M     is  singular  then  there  is  a  nonzero  vector  p     such  that 


M     p   -   0 
o  —  o 


(C-78) 
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M     is  symmetric  thus  it  is  also  true  that 


and 


This  implies  that 


pT   M     =  0  (C-79) 

—  o       o  ' 


pT  M     p      =0  (C-80) 

—  o       o  — o  v  ; 


tf 


T      A  ^^T      A  1 
o 
1 


0  =     J"      P       ABBT    ATp     dt  (C-81) 

t. 


However  equation  (C-81)  is  the  integral  of  the  square  of  the  length  of  a 
vector 

*f  2 

f       |   BT    AT  p  dt  =  0  (C-82) 

ll 

Thus  for  M     to  be  singular  it  is  necessary  and  sufficient  that  the  vec- 
tor integrand  of  (C-82)  be  zero. 

BT(t)    AT(t)  pQ    =  0  (C-83) 

for  all  t    <  t  <  tf. 

That  a  vector  p      exists  for  the  optimum  trajectory  such  that  equa- 
tion (C-83)  is  satisfied  may  be  shown  by  application  of  the  Pontryagin 
maximum  principle  for  unconstrained  control.     The  Pontryagin  maxi- 
mum principle  asserts  that  if  the  final  value,  S,  of  some  combination 
of  the  state  variables,  x,  are  to  be  a  maximum  or  minimum  with  re- 
spect to  the  control  variables,  y;  that  is, 

min  S  =  dT  xf  (C-84) 

then  the  desired  path  will  be  one  such  that  the  gradient  of  the  state 
velocity,  x,  in  the  vector  space  of  all  admissible  y, 

x  =  g  (x,  y,  t)  (C-85) 
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will  be  orthogonal  to  a  vector  p  which  satisfies  the  differential  equa- 
tions adjoint  to  x  =  g  (x,  y,  t)  and  such  that  the  final  value  of  p  is 

pf  =  -d  +  v{  (C-86) 

where  d.  =  0  for  components  of  the  state  vector  which  are  fixed  at  tf 
and  v.     -  0  for  components  of  the  state  vector  to  be  minimized  or  max- 
imized. 

Thus  to  satisfy  equation  (C-84),  it  is  necessary  and  sufficient  that 


where 


3  (g  (x,  y_,  t) ) 

3y 


x  =  g  (x,  y,  t) 


P 


9  (g(x,  y,  t)  ) 


3x 


Pf  =  -  d  +  v_ 


T 


(C-87) 

(C-85) 
(C-88) 
(C-86) 


A  proof  of  the  maximum  principle  is  given  in  Appendix  D. 

Application  of  equations  (C-86),  (C-87),  (C-88)  to  the  guidance 


roblem  in  this  thesis  reve 

als  that 

Pf  - 

-ffi 

'  •  {» 

T 

P =  - A1  p 

T 
P 

B  =  0 

(C-89) 


(C-90) 


(C-91) 

T 
where  A      and  B  are  the  matrices  of  partial  derivatives  in  (C-88)  and 

(C-87)  respectively. 

T 
Since  p  satisfies  equation  (C-90)  and  also   A     satisfies  (C-90),  a 

solution  of  p  is 

P=   ATPf  (C-92) 
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Then  substituting  into  equation  (C-91) 

p^  AB  -  0  (C-93) 

for  all  t,  <  t  <  t..     Thus  p.is  a  null  vector  of  M   .     Furthermore,  since 
1  i  —i  o 

t1  is  arbitrary 

MQ(t1)  pf  =0  (0<  tx  <  tf)  (C-94) 

To  show  that  M  is  not  singular  it  is  sufficient  to  show  that  every 
null  vector  of  M     satisfies  the  boundary  condition  on  p,  equation  (C-86). 
This  sufficiency  condition  is  easily  verified  by  noting  that  the  final 
element  of  pf  is  invariant  and  that  M  may  be  obtained  from  M     by  de- 
leting the  last  row  and  last  column  of  M   ,     In  general,  if  TV  is  inter - 

O  rp  rp 

preted  as  the  nonsingular  transformation  that  transforms  pf  into  p 
T  — i  — 

and  if  p  is  a  null  vector  of  B  because  of  the  component  p-  o,  then  delet- 
ing the  ith  row  of  TV  to  produce  A*  destroys  the  transformation  which 
carries  p. f  into  p. 

The  fact  that  only  the  last  row  of  A  need  be  deleted  is  a  con- 
sequence of  the  fact  that  only  the  last  component  of  the  final  state  vec- 
tor, mf,  has  been  maximized  in  obtaining  the  optimum  path.     Therefore, 
we  assert  that  if  every  null  vector  of  M     satisfies  equation  (C-89), 
the  M  obtained  in  equation  (C-77)  is  nonsingular. 

Let  us  prove  then  that  all  null  vectors  of  M     satisfy  the  boundary 
condition.     For  the  problem  in  this  thesis  B  is  a  seven  by  three  matrix 
of  rank  three  and  contains  zeros  in  the  first  three  rows.     By  virtue  of 
the  rank,  the  three  columns  of  B  are  linearly  independent.     Further, 
since  B  is  of  rank  three  there  are  four  linearly  independent  vectors 
which  satisfy  (C-91)  and  thus  are  null  vectors  of  B.     Three  of  the  four 
vectors  may  be  chosen  to  satisfy  (C-91)  by  virtue  of  the  three  rows  of 
zeros  in  B.     These  vectors  do  not  satisfy  (C-92)  and  (C-93)  and  there- 
fore are  not  null,  vectors  of  M   .     The  fourth  vector  must  be  a  null  vec- 

o 

tor  of  M     since  the  three  columns  of  B  and  the  four  null  vectors  of  B 
o 

form  a  basis  in  seven  dimensional  vector  space  and  six  of  the  vectors 
are  not  null  vectors  of  M   .     This  vector  is  unique  and  must  satisfy 
(C-89)  through  (C-93). 
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The  proof  is  thus  complete  since  all  vectors  may  be  written  as 
linear  combinations  of  the  vectors  in  the  basis. 
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APPENDIX  D 
PONTRYAGIN'S  MAXIMUM  PRINCIPLE 

D.  1  Introduction 

In  this  appendix  a  derivation  of  Pontryagin's  maximum  principle 
is  given  which  is  applicable  to  the  problem  of  the  thesis. 

Consider  that  the  cost  function  S  is  some  scalar  function  of  the 
final  state  variable,  xf.     It  is  always  possible  to  cast  a  problem  into 
this  form  by  redefining  variables.     The  control  vector  will  be  denoted 
by  a  generalized  vector  u,  which  is  bounded  by  certain  constraints.     In 
addition,  the  solution  must  satisfy  boundary  conditions  on  some  of  the 
state  variables.     The  desired  solution  is  the  optimal  control,  u  ,  which 
will  maximize  the  cost  and  satisfy  all  constraining  conditions. 

D„  2  The  Maximum  Principle 

Assume  that  the  state  variables  may  be  written  in  the  form 

X  =  £(X,  u)  (D-l) 

and  that  the  solution  is  to  satisfy  the  final  boundary  conditions 

£(xf)  =  £  (D-2) 

The  cost  function  may  be  written  as 

S  -  dT  xf  +  vT  £  (D-3) 

where  d  is  a  known  vector  and  v  is  an  unknown  constant  vector. 

Then  small  changes  in  the  cost  due  to  changes  in  state  are  given  by 

6S  =  d     6xf+/    6xf  (D-4) 

3    X    n 


where 


d  (j)  . 

-  1  or  0  (D-5) 

3    X.    n 

if 
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Equation  (D-5)  is  1  for  fixed  boundaries;  0  for  unspecified  boundaries. 
The  variables  adjoint  to  (D-l)  are  given  by 

p  =   -  A'p  =   -  I  j       p  (D-6) 

P  (tf)    =  pf  (D-7) 

The  vector  p  is  frequently  called  the  "costate". 

Assume  that  an  optimal  control,  u  ,  is  known.     Then  for  small 
changes  6u,  the  state  must  change  by  6x.     Then 

6x   =  g  (x°  +  6x,  u°  +  6u)   -  g  (x°,  u°)  (D-8) 


Consider  the  term 


—  (pT6x)  -  pT6x  +  pT8x  (D-9) 


dt 

Applying  (D-8)  one  obtains 

"To          d      /Tr>  Tr      /O.c  o,rx         ,   o       oxn 

p     ox  = (p     ox)  -  p     [g  (x     +  ox,  u     +  ou)  -  g  (x  ,  u    )  ] 

dt  (D-10) 

Integrating  (D-10)  yields 

J    p     6x  dt  =  p     6x  ~     /    P      [g  (x°  +  5x3  u°  +  6u)  -  g  (x°,  u°)]  dt 

0  0        0 

(D-ll) 

If  the  initial  state  boundary  is  fixed,  6x      =0  and  the  first  right  hand 
term  is 

6S  =  pj   6xf  (D-12) 

provided  pf  is  chosen  such  that 

T    T    T  d  $- 

Pf  =  d  +  "      (D-13) 

3  x„ 

The  right  side  integrand  may  be  expanded  in  a  Taylor  series  around  the 
optimal  trajectory  if  the  variations  due  to  control  and  due  to  state  are 
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separated.     That  is 

g  (x     +  6x,  u     +  6u)  -  g  (x   ,  u    )  =  g  (x   ,  u     +  6u)   -  g  (x   ,  u   ) 

+  [g  (x  ,  u     +  6u)  ]  6x  +  remainder 

9-  (D-14) 

If  the  adjoint  relationship,  (D-6),  is  used  on  the  left  hand  integrand 
of  (D-ll),  the  term  may  be  rewritten  as 

/    pT  6x  dt  =  -    /    pT  —[g  (x°,  u°)  ]  6x  dt  (D-15) 

0  0  3x 

Then  substituting  (D-12),  (D-14)  and  (D-15)  into  (D-ll),  the  result  is 


*f 


&  S  =    /     p1  [g  (x°,  u°  +  6u)  -  g  (x°,  u°)  ]  dt 
0 


tf 

T 

(ax  ,      J 


J     pT\—  [g  (x°,  u°  +  6u)  -  g  (x°,  u°)]  )    6x  dt     (D-16) 


Consider  the  first  integrand  in  (D-16),     If  S  is  a  maximum,  then  any 
allowable  change  6u  must  cause  6  S  to  be  negative,     Therefore  for  all 
0  <  t  <  tf 

pT  [g  (x°,  u°  +  6u)  -  g  (x°    u°)  ]  <    0  (D-17) 

In  considering  the  second  integrand  it  is  argued  that  admissible 
control  changes  6u  ,  can  produce  only  a  small  variation  in  the  deriva- 
tive term  such  that  its  product  with  6x  is  second  order  or  higher.     Con- 
sequently the  second  integral  and  all  higher  order  remainder  terms  of 
(D-14)  are  neglected. 

Therefore  (D-17)  represents  the  sufficient  condition  for  an  optimum. 
A  necessary  condition  is 

pT  [g  (x°    u°  +  6u)  -  g  (x°,  u°)  ]  <    0  (D-18) 

The  implication  of  (D-18)  is  that  for  all  points  on  the  optimal  tra- 

T 

jectory  and  for  all  admissible  functions  u,  p    g  must  be  a  maximum. 

Replacing  g_  with  x,  notice  that  the  scalar  product  of  p  and  the  state 
velocity  must  be  a  maximum 
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If  the  control  is  unconstrained,  then  a  necessary  condition  that 

T 
p    |is  a  maximum  is  that 

a       t  t   9g- 

(P     g)  =  0  =  P       (D-19) 


3  u  3  u 

From  previous  work 


B  =   (D-20) 

3u 

Thus  a  necessary  condition  for  a  maximum  is 

pTB  =  O  (D-21) 

It  is  possible  to  derive  the  Hamiltonian  formulation  used  in  Chapter 
III  by  introducing  a  slightly  different  cost  function  S  ,    In  place  of  (D-3) 
let 

.T 


S'  =    /    h  (u)    dt  +  v     £  (D-22) 

0 


where  h(u)  is  some  function  of  the  control.     Then 

tf 

/  _      _ 

0  3u  3x 


6S'  :     J       —    6udt+^i  6x  (D-23) 


■f 

where  again 

30. 

—    =   1  or  0  (D-5) 

3  x  .  p 
if 

Equations  (D-6)  through  (D-ll)  hold  and  (D-23)  may  be  written  as 

tf 

6  S1  =    /    [h  (u°  +  6u)  -  h  (u°)  ]  dt  +  pT  &xf  (D-24) 

0 


provided  that  pf  is  chosen  such  that 


rp  rp     3^ 

pj    =  v1  (D-25) 

3x» 
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If  D-24)  is  substituted  into  (D-ll)  along  with  (D-14)  and  (D-15)  ,  the 
result  is 

6  S'  =    /       pT  [£  (x0,  u°  +  6u)  -g  (x0,  u°)  ]  +  [h  (u°  +  6u)  -h  (u°)]  \  dt 

o      I  J 

L  u-    u  a       *  (D-26) 

+  higher  order  terms  v  ' 

Equation  (D-26)  is  analagous  to  (D-16),     A  sufficient  condition  that  S' 
be  a  maximum  is  that  for  all  admissible  control  changes,   6u 

pT  [g  (x°,  u°  +  6u)  -  g  (x°5  u°)  ]  +  [h  (u°  +  Su)  -h  (u°)  ]  <  0    (D-27) 


The  implication  of  (D-27)  is  that  for  all  points  on  the  optimal  trajectory 
and  for  all  admissible  functions  u°,  (p    g  +  h(u)  )    must  be  a  maximum. 
Using  the  VSI  vehicle  as  an  example 

2 


Then  using  (D-l)  define 


h(a)    -  —  (D-28) 


a2  T 

H  -  —  +  p1  x  (D-29) 


For  the  linearized  guidance  problem  of  Chapter  III,  the  vector  ^ 
was  interpreted  as  the  velocity  of  the  state  error  at  the  final  time. 
Since  i[  is  a  state  velocity  and  pf  is  the  final  value  of  the  costate,  using 
(D-25) 

a2  T   • 

H  =  —  +  V     |  (D-30) 


where  for  fixed  boundary  conditions 


34 


(D-31 


Equation  (D-30)  is  a  linearized  version  of  the  general  Hamiltonian 
formulation    but  with  a  nonlinear  cost. 

D.  3  Optimization  Criterion 

It  was  briefly  mentioned  in  Chapter  V  that  the  optimality  criterion 
for  Pontryagin's  maximum  principle  is  more  workable  than  the  calculus 
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of  variations  criterion.     This  remark  deserves  further  discussion. 
From  the  previous  section  it  may  be  noted  that  the  criterion  for  optim- 
ally may  be  stated:   "For  admissible  control  variations,  the  change  in 
the  cost  function  must  be  zero  or  negative  (positive)  if  the  cost  is  max- 
imum (minimum).     Admissible  control  variations  are  those  which 
violate  neither  the  boundary  conditions  nor  any  control  constraint.  " 

In  contrast  to  this,  the  calculus  of  variations  criterion  for  optim- 
ality  is  derived  from  the  integral  of  the  functional  F  used  in  Chapter  V. 
Referring  to  equation  (5-18),  for  example,  the  optimality  criterion  may 
be  stated:   "For  admissible  control  variations  and  variations  in  state, 
the  boundary  conditions  must  be  satisfied  and  the  cost  must  not  change 
to  first  order.  " 

Since  the  integral  of  the  functional  equals  the  cost,  the  form  of 
equation  (5-18),  rewritten  here  as  (D-32), 

r       itf     tf 

S/Fdt=0=  +  /      [(         )Srn-(        )  6y_  + 

JO       0 

(  )  6m  +  (         )  6a       dt  (D-32) 

gives  the  impression  that  the  criterion  must  hold  for  all  arbitrary 
variations  in  state  and  control.     The  implication  is  not  true  for  the 
problem  of  this  thesis  and  is  misleading  at  best  because  (D-32)  does 
not  represent  the  first  variation  of  the  cost.     It  will  now  be  shown 
that  although  (D-32)  is  satisfied  for  all  variations  in  state  and  control, 
the  optimality  criterion  can  not  be  satisfied  for  all  arbitrary  state 
variations. 

Let  the  guidance  equation  (3    27)  be  partitioned  such  that  it  con- 
sists of  two  equations 

tf 


6  r 


(D-33) 


tf 


6m     =    A  *   6s  t  +    /     TV  ^  B  6a  dt  (D-34) 
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Using  the  expansion  of  the  matrix  _/\_ ,  (D-33)  and  (D-34)  are 


8  r 

6v 


An     A12 

A21     A22 


*f  ' 


/  6-  1     '+      /     /V  B  6a    dt 
\6v/  t 


(D-35) 


1  2 

6mf=     A336mt-    /      A33   (?L)      aT  Sa  dt 

t  P 


(D-36) 


2        2 
Using  the  solution   Aoo  =  mf/m     from  equation  (5-25),  (D-36)    maybe 


33 


reduced  further 


m 


8m 


f  2 

mt 


2  tf 

m  1 

*          f  r     T  a   a, 

J       a      8a  dt 


t 


(D-37) 


P      t 


For  the  criterion  to  be  satisfied  it  is  necessary  that  the  final  state  varia- 
tion 8s  n  vanish.  First  allow<  —  >   to  take  on  a  nonzero  value  but  let  8m, 
~f  l&v;  t 

equal  zero.     Then  (D-35)  and  (D-37)  become 


O  =  £  +    /    A'  B  6a  dt 

t 


m 


2        t 


0     =     0- 


<-f 


f 


a      8a  dt 


(D-38) 


(D-39) 


P      t 


From  Chapter  III  it  is  known  that  all  propellant  optimal  control  varia- 
tions which  satisfy  the  boundary  conditions  are 


6a  =  BTA    *T  M"1  (22  -  |)  -a 

Substituting  this  solution  into  (D-39)  and  expanding  yields 

2 

7]TM_1?7-  2J  -^M"1  £ 


P      L 


However,  along  an  optimal  trajectory 


7]T  M'1  r]  -  2J  =  0 


(D-40) 


(D-41) 


(D-42) 
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Thus  to  satisfy  (D-39)  it  is  necessary  that 

0  -  ^M"1  |  (D-43) 

(D-43)  does  not  hold  in  general.     It  holds  only  when  the  error  |  is  such 
that  a  control  change  6a  orthogonal  to  a  can  satisfy  boundary  conditions. 

Now  consider  a  variation  in  6m  „  but  no  variation  in  position  or 
velocity.     Then  it  is  required  that 

tf 

O  =     /    A*  B  6a  dt  (D-44) 

t 

2  2  t 

m  f         m  „ 

0  =  — *-     -  — L      J     aT  6a  (D-45) 

mt  p       0 

A  general  variation  6a  may  always  be  written  as 

6a  =  BT    A*T  I  -  a  (D-46) 

Inserting  (D-46)  into  (D-45)  and  (D-44)  and  solving,  one  obtains  from 
(D-44) 

n  =  M_1  T)  (D-47) 


and  from  (D-45) 


6m 

p  ~  =  77     w  -  2J  (D-48) 

mt 


Eliminating  n 


6mt  T      -1 

p j-     =  ^     M      2Z  -  2J  =  0  (D-49) 

m  „ 

Therefore  there  are  no  variations  6m     which  permit  the  boundary  con- 
ditions to  be  satisfied  and  which  do  not  change  the  cost. 

It  is  concluded  that  the  calculus  of  variation  criterion  and  formula- 
tion can  be  misleading  if  not  treated  with  much  care.     The  Pontryagin 
principle  avoids  the  difficulty  by  making  no  explicit  requirements  on 
state  variations. 
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APPENDIX  E 
STEEPEST  ASCENT  FORMULATION  FOR  GUIDANCE 

E.  1  Introduction 

In  this  appendix  the  guidance  problem  will  be  formulated  as  it  would 
appear  in  a  "steepest  ascent"  trajectory  computation  procedure,     The 
steepest  ascent  method  is  characterized  by  a  rigid  control  over  the 
step  size  between  iterations.     This  step  control  is  to  prevent  the  pro- 
cedure from  violating  linearity  assumptions  as  it  approaches  the  opti- 
mum.    Sophisticated  techniques  have  been  devised  for  automatic  selec- 
tion  of  step  size.  However,  at  some  point  the  investigator  must  use 

his  experience  and  judgement  to  select  a  constant  or  matrix  of  con- 
stants to  insert  in  the  selection  procedure  for  step  size,      The  author 
is  convinced  that  the  additional  complexity  of  rigid  step  control  cannot 
always  be  justified. 

E.  2  Steepest  Ascent  Solution 

For  simplicity  in  illustrating  the  steepest  ascent  method,  consider 
only  the  trajectory  problem  for  the  unrestricted  VSI  mode  of  control. 
In  this  problem  it  is  desired  that  the  optimal  acceleration  integral  be 
a  minimum. 

tf    ,       ,2 
1        o 

J°  =     J       '-    '      dt  (E-l) 

0  2 

It  is  further  desired  that  the  change  in  the  acceleration  program  from 
one  iteration  to  the  next  be  such  that 

tf  2 

k2  =     /     |  6a  |      dt  (E-2) 

0 

2 
where  k     is  a  constant  to  be  selected  at  the  discretion  of  the  investigator. 

In  order  to  satisfy,  to  first  order,  the  terminal  constraint     jj,   6a 
must  satisfy  the  relation 
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-  £  =  /     X'b  6a  dt 

0 

where  all  variables  are  defined  exactly  as  in  Chapter  III. 
Define  the  new  cost  function  J'  such  that 


(E-3) 


T 
tf      O x     O 

1  a      a 


0  2 


dt  + 


tf  tf 

— (k2   -   /    6iT  6a.  dtV  12    (  i+  /      ^*B  6*  dt) 

(E-4) 

where  it.   and  7r  „  are  constant  Lagrange  multipliers,  a  scalar  and  a 
vector  respectively. 

For  arbitrary  variations  in  6a,,  the  first  variation  of  J'  must  vanish. 
Thus  for 


o 


a     -  a  +  6a 
where  a  is  a  nonoptimal  program, 


0  = 


(a  +  6a)T  =  tt1  6aT  +  j^   A"~B 


(E-5) 


6  (6a)  (E-6) 


5  (6a)  /  0 


or 


6a  = 


T       *T 
-B1    A       £2  "  - 

(1   "  JrJ 


(E-7) 


The  multipliers  may  be  determined  by  satisfying  the  constraint  equa- 
tions.    First  solve  for  ir  9  through  equation  (E-3). 


jr     =  M"i[(l-7r1)  |  -  7?  ] 


(E-8) 


where 


T       .*T 


M  =      f    A!itBB        A.       dt 
0^ 


(E-9) 


T7  =     /      A*  Ba  dt 

0 


(E-10) 
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and 


.  T 

3f     J- 


"I 


I         " 


6  a  =  -  B      A"      M 
Using  (E-2)  to  solve  for  one  obtains 


n 


(i-^)J     (i-^) 


(l-7T1) 


(i-^r  = 


2 J  -  ?7TM~% 
k2  -  |TM_1| 


or 


(1-^) 


k2  -  |TM_1| 

T      ^1 
2  J  -  r\     M      77 


(E-ll) 


(E-12) 


(E-13) 


where  J  is  the  acceleration  integral  for  the  nonoptimum  program. 


The  factor 


(i^r) 


is  the  step  size  control  which  must  be  deter- 


mined by  judicious  selection  of  k   .     Observe,  however,  the  result  of 

T 
using  (E-ll)  to  evaluate  the  integral  of  a     Sa.     This  term  is  the  cross 

product  in  the  expansion  of     |  a    |       by  equation  (E-5). 


r      T                         T      - 1              1  T      - 1 

J    a     6adt=-?)     M      { —  (2J  -rf  M      77) 

0  1-7T. 


(E-14) 


1 


o 


As  the  program  a  approaches  a  ,  the  difference,   6a,  must  vanish. 

Likewise  the  error  £_  must  vanish.     From  (E-14)  either  the  step  size 

T      -1 
must  approach  zero  or  the  term  (2J  -  77     M      77  )  must  approach  zero  or 

both.     It  is  easy  to  show  that  both  must  approach  zero.     Consider 

a    =    BT  A*TM-1?7    -a  (E-15) 

Then 

tf 
/    aTa  dt  =  2 J 


77     M      77 


(E-16) 
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Clearly 

f      aTa  dt  >  0  (E-17) 

o 

with  equality  holding  only  when 

a  =  BT  A*TM"%  (E-18) 

Equation  (E-18)  holds  only  along  the  optimal  trajectory.     Therefore  from 

/     1     \  2         T      - 1 

(E-13)  it  is  clear  that  for  i  ■=— — J  to  remain  finite,  (k     -  £     M       |)  must 

T      - 1  II 

approach  zero  at  least  as  fast  as  (2J  -  r)     M      77).     Since   ||  |  approaches 

2  ~~ 

zero  this  implies  that  k     must  be  reduced  to  zero  as  the  optimal  path  is 

approached.     One  concludes  that  "steepest  ascent",  as  formulated  using 

the  usual  technique       ,  inherently  converges  slowly. 

The  approach  in  this  thesis  is  to  assume  that  as   K  I  becomes  small 

the  step  size  will  automatically  become  small  without  the  arbitrary  con- 

2 
straint  imposed  by  k   .     Thus  ir  ^  is  set  equal  to  zero  and  the  procedure 

allowed  to  converge  as  rapidly  as  possible.     The  results  appear  to 

justify  this  procedure. 
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APPENDIX  F 
COMPUTATIONAL  COORDINATES 

F.  1  Introduction 

In  this  appendix  the  computational  coordinate  system  is  derived 
from  the  ephemeris  data  of  the  launch  planet  and  target  planet.     In  the 
chapters  describing  the  guidance  equation  and  its  uses,  generalized 
vector  notation  is  used  and  the  problem  of  coordinates  does  not  arise. 
The  motion  is  assumed  to  be  described  in  a  nonrotating  frame. 

For  computations,  however,  it  is  necessary  to  be  more  specific. 
The  computational  frame  used  in  computer  tests  is  a  heliocentric  frame 
defined  by  the  transfer  plane  and  the  initial  point.      The  transfer  plane 
is  the  plane  which  contains  the  initial  point,  the  final  point  and  the  sun. 
The  x  axis  of  the  system  passes  through  the  initial  point,  the  z  axis  is 
northerly  and  the  y  axis  completes  the  triad. 

The  objective  is  to  describe  the  transfer  plane  and  its  coordinate 
system  in  terms  of  ephemeris  data.     This  data  may  be  given  in  either 
ecliptic  or  equatorial  coordinates.      The  description  of  these  systems 
and  the  transformation  between  them  is  available  in  most  basic  celes- 
tial mechanics  texts  and  will  not  be  reproduced  here. 

F.  2  Computational  Coordinates 

Designate  the  unit  vectors  in  any  system  as  i,  _j,  and  k  and  attach 
subscripts  to  denote  the  system.     For  the  computational  system  use 
subscript  c.     Then,  if  the  launch  and  target  points  are  denoted  by  sub- 
script L  and  T  respectively, 

-T 

rL 

r  T    X  r  _, 

k      =    — —  (F-2) 

|rLXrTl 

J      =  k     X  i  (F-3) 

—  c      —  c     —  c 
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The  transfer  angle,   9 ,  in  the  transfer  plane  is  determined  by 


sin  9   - 
c 

rLXrT 
rLrT 

k      sin  9 
—  c 

=  sin  9 

(F-4) 


(F-5) 


The  coordinates  of  the  launch  (initial)  point  are 


XL 


y 


L 


'L 


L 
0 

0 


(F-6) 


The  coordinates  of  the  target  point  are 

/       \ 


r^  cos  9 


*T 


rT  si 


n  9  \ 


(F-7) 


Two  additional  parameters  of  interest  are  the  inclination  of  the 
transfer  plane  to  the  planes  of  the  launch  planet  and  the  target  planet. 
If  the  northerly  normals  to  these  planes  are  denoted  by  kT    and  kT 
respectively,  and  the  angles  by  a      and  a      respectively,  then 


cos  aL  =    kL 


(F-8) 


cos  a^  =    k  „  '    k 
r       —  1       — c 


(F-9) 


where 


isL 


-L  X-(L  +  90°) 
-L  X-(L  +  90°) 


(F-10) 


_  ry 


£T  Xr^  +  9qo^ 

r       X  r 

LT      -(T  +  90°) 


(F-ll) 
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From  the  preceding  equations,  numerical  values  of  the  coordinate 
transformations  (i.  e.    the  orthogonal  transformation  matrix)  between  an 
ephemeris  tabulation  and  the  computational  coordinates  can  be  com- 
puted when  the  launch  date  and  transfer  time  are  specified. 
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APPENDIX  G 
COMPUTER  PROGRAM 

G.  1  Introduction 

The  FORTRAN  program  used  to  test  the  guidance  theory  is  repro- 
duced in  this  appendix.     The  input  data  format  is  given  as  well  as  the 
units  for  input  data. 

There  are  a  few  areas  in  the  program  which  can  be  made  more 
efficient,  since  testing  is  completed,  however  it  operates  quite  satis- 
factorily as  presently  written.     The  author  is  indebted  to  Mr.    Krupp 
for  his  magnificent  efforts  in  writing  this  program. 

G.  2  Input  Format  and  Units 


Nine  data  cards  are  required  for  each  problem  to  be  computed. 
The  format  for  each  is  3F20.  9.     The  following  sequence  and  units  are 
required: 

1.  )  Initial  position  (A.  u.  ) 

2.  )  Initial  velocity  (A.  u.  /day) 

3.  )  Target  position  (A.  u.  ) 

4.  )  Target  velocity  (A.  u.  /day) 

5.  )  and  6.  )      Estimated  v_  vector.     For  VSI  control,  both  cards  are 

null  vectors.     For  CSI  control,  card  number  5  is  the  null 

vector  and  card  6  should  contain  a  number  of  the  order  of 

2 
maximum  initial  acceleration.     The  units  are  (A.  u.  /day   ). 

7.  )       a.     Number  of  iterations  desired 

-4 

b.      Maximum  initial  acceleration  in  units  of  10      g   .     (Example, 

1.2).     For  VSI  problems  use  any  large  number. 


c.     Gravitational  constant  of  the  central  body  in  units  of 

3  2 

(A.  u)   /day   .     For  the  sun  this  vali 

)       a.     Maximum  time  step  desired  (days) 


(A.  u)3/day2.     For  the  sun  this  value  is  0.  000295912. 
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b.  Earth  gravity  (A.  u.  /day   ).     This  value  is  0.  494840. 

c.  Flight  time  (days)  . 

2  3 

9.  )       a.     Exhaust  power  of  engine  per  unit  mass  (A.  u.  )   /day  . 

/  4 

The  conversion  is:  divide  power  in  kw/kg  by  3.  4653  X  10    „ 

b„     Switching  number,      (use  1.  0), 

The  routine  will  accept  as  many  problems  of  one  type  as  desired.     CSI 
and  VSI  problems  cannot  be  run  together.     Submit  a  complete  set  of 
data  cards  for  each  problem. 

G.  3  The  FORTRAN  Program 

Except  for  subroutine  DERV,  the  programs  for  CSI  and  VSI  control 
are  identical,,     There  is  a  subroutine  DERV  for  each  mode  of  engine 
control.     Select  the  routine  appropriate  to  the  problem  and  omit  the 
other.     Selection  cannot  be  made  automatically  with  the  current  pro 
gram. 

The  FORTRAN  source  program  is  reproduced  on  the  following 
pages.     The  generalized  flow  chart  is  presented  in  Figure  G-a.     Figure 
G-b  is  the  storage  map  for  computer  variables. 
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Fig.  G-b.  Computer  storage  map. 
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LABEL 
LIST 

M  A  I  N 

D  I  MENS  I ON  W ( 1 2  »  2  0),  T ( 7 ,  8  » 

READ  1»  ( (W( I »  J) »  J=15»  20) 

INPUT  DATA,  (SEE  FORMAT) 
LAUNCH  (POSITION,  VELOCITY) 
TARGET  (POSITION,  VELOCITY 
ESTIMATED  NU  VECTOR 

NUMBER  OF  ITERATIONS.  MAX  ACCELERATION, 
MAX  TIME  STEP,  EARTH  GRAVITY,  FLIGHT  TIME 
POWER,  SWITCHING  NUMBER 
FORMAT ( 3F2  0.9 ) 
NIT=W(7,  15) 

W(S,  20)=W(7,  19  1/10000.0 
W(6»  17)=W(7,  16-)*W(8»  20) 

TRANSFER  OF  INPUT  QUANTITIES 
W(8»  16)=2.0*W{8»  17)/W(6,  17) 

COMPUTATION  OF  EXHAUST  VELOCITY 
HP  5    L=l-  NIT 

LIMIT  (NIT ) 

T  ) 

T) 

I  )  ,  1=15,  20) 


7  ,  2 )  ,  W ( 8  »  17),  W ( 8 ,  15) 


GRAVITATION  CONSTANT 


FO  FINAL  STORAGE 


I TERATION 
CALL  CORRECT*', 
CALL  LAMTAB(W, 
PUNCH  2,  (W(5, 
FORMAT ( 3E20.8 ) 
GO  TO  3 
END 


27 
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10 


LAMTAB 

LABEL 

LIST 

SUBROUTINE  LAMTAB (W,  T) 

SUBROUTINE  COMPUTES, 
DIMENSION  W( 12,  20  )  ,  T  (7  , 
L=V(8,  19) +.01 

L  IS  STEP  NUMBER 
DO  2^  LL=1  ,  L 
A  =  T ( 7 ,  1 .  L  L) * T ( 7  ,  1  , 
IF( 5*(  ( LL-1  ) /5 1-LL+]  ) 

PAGE  SPACING  OF 
PRINT  3 


PRINTS  BACKWARD  INTEGRATED  LAMBDA 

8,  60  )  ,  Q( 10»  10) 


MATRIX 


L+l  ) 
6>  5, 

LAMBDA 


6 
PRINT 


OUT 


DO 
no 
0(  I 
DO 
0(  I 


10  I: 
10  U  = 
»  J)  = 

10  K  = 


J)=Q(I,  J)+T(I,  K.+  1,  LL)*T(K»  U+l,  L+l) 
INVERTED  FINAL  VALUE  OF  FORWARD  INTEGRATED 
MULTIPLIED  INTO  LAMBDA  MATRIX  AT  EACH  STEP 

FORMAT(104X>  F16.8) 

FORMAT ( 1H12X4HTIME50X13HLAMBDA  MATRIX ) 

FOPMATdH  ,F7.2,  7F16.8) 

FORMATI8X,  7F16.8) 

PRINT  2 

PRINT  2,  T(l,  1,  LL),  (Q(I,  1).  1=1,  7) 

PRINT  1  ,  ((0(1,  J),  1=1,  7),  U  =  2 ,  7) 

PRINT  4  ,  A 

RETURN 

END 


LAMbDA  MATRIX 


30 
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COQPEC 

LABEL 

LI^T 

SUBROUTINE  CORREC(W,  T) 

SUBROUTINE  COMPUTES  NEW 
AMD  PRINTS  FINAL  VALUES 

DIMENSION  W(12.  20),  Q(10»  10) 

CALL  INTEG(W,  T) 

DO  10  1=1,  6 

V  (  R 

W(  4 


INITIAL  VALUES  FOR  NEXT  ITERATION 
OF  INTEGRATED  QUANTITIES 
T (7,  8,  50) 


AT  TARGET  (TARGET  XI ) 


( 3,  U  +  14) ) 


1 3)  =2.0 

1  +  14) =W<  1  »  I  )-W( 3,  1+14) 
COMPUTATION  OF  MISS  VECTOR 
W  (  2  »  1  +  14  1=0.0 
DO  10  J  =  1 ,  ft 
W(2,  1  +  14)  =W(  2,  I+141+WU+5,  J)*(W(1»  J) -I 

TRANSFORMS  TARGET  XI  TO  EQUIVALENT  LAUNCH  ERROR  (INITIAL  XI) 
10  0(1,  U ) =W( 1+5 ,  J+7 ) 
CALL  INVERT(6,  0) 

INVERSION  OF  M  STAR  MATRIX 
DO  20  1=1,  6 
DO  2  0  J=l ,  6 
2 'J  W(5»  1  +  14)  =W(  5,  1  +  14) -0(1.  J)*W(2t  J+14) 
COMPUTATION  OF  CORRECTED  NU  VECTOR 
3  FORMAT ( 1H15X13HINITIAL  STATE5X 1 1HF I NAL  ST ATE6X 1 2H TARGET  STATE9X 

C  3HETA1 1X9HTARGET  X  I  8 X 1  OH  I N I T  I  AL  XI7X10HNU  X  10000) 
2  FORMAT (1H  ,7F17.<?) 

1  FORMAT ( 1H0/1H09X1HJ5  3X12HMSTAR  MATR I X/1H0EI7. 7 »6F1 7, 
AJ=W(5»  3)/2.0 
P  R  I  M  T  3 
P  R  T ,'  1  T  2 

DO  2  5  1=15-  20 
XNU=W( 5 ,  I  (-10000. 0 
2  5  PRINT  2  ,  W (  1  ,  I  )  ,  W ( 1  ,  I  - 1 4 )  , 
PRINT  1  ,  A  J  ,  (  (  ,v  (  I  ,  U  )  ,  I 
DO  3  0  1=1,  7 
no  3  0  J=l  ,  7 
30  Q(  I ,  J) =W( J  +  5,  I  ) 
CALL  T  N  V  E  R  T ( 7  ,  Q ) 
L  =  W  (  R  ,  1  a  )  +  1  .  0  1 
DO  40  1=1,  7 
DO  40  J=2,  8 
4"  T(  I  ,  J,  L)=Q(  If  J-l ) 
T(7,  1,  L)=1.0/W(5, 


6/(  loX6F17.6)  ) 


/  (  3  ,  I  )  ,  W (  3  ,  I  -  1 4  )  ,  W  (  4  ,  I  )  , 
6,  11),  U  =  8  ,  13 ) 


(2,1  )  ,XNU 


7) 

NVERSION  OF  FINAL  VALUE  Or  FORWARD  INTEGRATcD  LAML-DA  MATRIX 


RETURN 
FWD 
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IMTEG 
LA3EL 
LIST 
SUBROUTINE 


INTEGIVi,  T) 


C  SUBROUTINE 

C  AND  PRINTS 

DIMENSION  W(  12  , 
DO  10  1=1,  168 
D  (  I  ,  1  )  =  D  .  n 

i  o  w  ( i ,  i  j  =  o .  a 

DO  2  3  1=1,  6 
W  (  I  +  5  ,  I  )  =  1  .  0 

I  )  =  W  (  1 »  I  + 
I ) =W( 5. 
7) =1.0 
7 ) = j.C 
7)  =1.0 
7  )  =  1  .  0 


CONTROLS  INITIAL  VALUZo  Fuk 
OUT  TRAJECTORY  VARIABLES  AT 
20  )  ,  D( 12  »  20)?  T ( 7,  8,  50) 


INTEGRATION 
EACH  STEP 


(  1  , 
(2  > 
(  1  . 
(2, 
(  5  , 
(  12 


14) 
1  +  14) 


■•/(  7 
>/  (  7 
•M  3 


1  9  )  =  0  .  0 

1  5  )  =  1  .  0 

18  )=: 

IN  I TI AL 
DERV( W» 


FOR  INTEGRATION 


32 


CONDI  TIO-MS 
CALL  DERV( W »  D ) 
I F ( V ( 7  ,  1 6  )  -W ( 6  »  17))  30,  3  0  »  35 

TEST  FOR  THRUST  LIMITATION  SWITCH  POINT 
'<'  (  7  ,  1  5  )  =0  .  0 
CALL  DERV(W,  D) 
D  E  L  = W { 7 ,  16) 

LIMITATION  OF  INTEGRATION  oTEP  SIZE 
FORMAT ( lhlliX4HTIME8X12nACCELERATIQN7X9HACCEL  MAG8X8HPOS1 TION9X 
C  8HVELOCITY11X4HMAGS13X5HGAMMA/26X9HX  10000/G8X9HX  10000/G) 
FORMAT (1H  ,F17.4,  6F17.7) 
FOPMATI1H  ,17X,  F17.7,  17X,  2F17.7) 
DO  1-)     i  =  i,  5  0 


4  5 
46 


6 

7 'J 

8  0 


A  1  =W(  6,  18)  /v.MR,  20) 

A  2  =  W ( n  ,  19)/ W (  8  ,  20) 

A3=l\'(6»  20)/W(8»  20) 

AM=D(5»  2)/'.\'(8»  20) 

GAM=W(7>  16)/W(6»  17) 

W(8»  1 9  )  =  I 

DO  4  0  J=],  8 

DO  40  K=l,  7 

T  (K  ,  J  ,  I  )  =  '•'  (  J  +  4,  K  ) 

I  F  (  1J#  (  {  I  - 1.)  / 1  0  )  - 1  +1  )  46  ,  4  5,  46 

PAGE  SPACING  FOR  TRAJECTORY  DATA 
p  ?  1 1\>  T  3 
p  P  T  M  T  2 

PRINT  2  ,  W ( 5,  1  )  ,  Al  ,AM,W ( 1 , 1  )  »W( 1  ,  4  )  ,  W (  1  ,  7  )  ,GAM 
PRINT  1  ,  A  ?  ,  W  (  1  ,  2  )  ,  V,'  (  1  ,  5  )  ,  A  3  ,  W  (  1  ,  3  )  ,  W  (1,6) 
TF(DKL+W(5,  1)-W(7,  20))  60,  60,  50 

TEST  FOR  END  OF  FLIGHT  TIME 
DEL=W(7,  20)-W(5i  1) 
IF(DFL/W(7,  18)-. 0001)  80,  80,  70 

TEST  FOR  FINAL  STEP  LESS  THmN  .0001 
CALL  S T  E  P ( W ,  D  »  DEL) 
R  F  T  U  R  M 
END 
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10 


1  o 

2  0 

25 

3  0 

35 


45 

5  J 
55 

6  0 

70 

75 


90 


INTEGRATION 
DK240).  W2(240)»  02(240).  W3(240),  03(240) 
FLAG 


STEP 

LABEL 

LI^T 

SUBROUTINE  STEP(W1»  Dl,  DEL) 

SUBROUTINE  SELECTS  STEP  SIZE  ON  dASIS  OF 

ERROR  AND  THRUST  LIMITATION  SWITCH  POINT 
DIMENSION  ,11  (  240  )  » 
FLAG=W1 ( 175 ) 

DEF  INITION  OF 
DC  10  1=1  j  240 
D  2  (  I  )  =  0 1  (  I  ) 
D  3  (  I  )  =  D 1  (  I  ) 
W  2  (  I  )  =  W  1  (  I  ) 
W3 (  I  )=W1(  I  ) 

DEF  INIT ION 
CALL  RUNKUT(W2» 

TRIAL  INTEG 
I  F  (FLAG- 'J 2  (  186  )  ) 


OF  VARIABLES 
D2,  DcL) 
RATION 

20,  15,  15 


FOR  TRIAL  INTEGRATION 


FOR  THRUST  LIMIT  SWITCH  OFF  POINT  DURING  INTEGRATION 


DURING  INTEGRATION 


1 ( 187)-W2(187) ) ) 

AT  SWITCH  POINT 


C^LL 
CALL 


TEST 
FLAG =0.0 
GO  TO  30 
IF (FLAG-MIN1F (W2 ( 186 )-1.0,  .5))  25,  40,  40 

TEST  FOR  THRUST  LIMIT  SWITCH  ON  POINT 
FLAG=] .0 
DFL=DFL#A3SF( (Wl ( 187  )-Wl  (  198  )  )  /  ( 

SELECTION  OF  STEP  BEGINNING 
DO  3  5  1=1,  ?40 
D  ?  (  I  )  =  D  1  (  I  ) 
W2 ( i  )=W1 (  I  ) 
CALL  R  U  N  <  U  T ( W  2  »  D 2 ,  DEL) 

INTEGRATION  OVER  FULL  STEP 

RUNKUTCW3*  D3,  DEL/ 2.0) 

RUNKUT(W3,  D3,  DEL/2.0) 

T'»'0  INTEGRATIONS  OF  HALF  STEP  EACH 
TEST  =  ABSF( M2 (  174  )  /W3 ( 174 ) -1 . 0 )/. 000045+ . 000 1 
IF(TEST-2.0)  55,  55,  50 

DIFFERENCE  TEST  OF  FULL  STEP  AND  HALF  STEP  INTEGRATION 
DEL=DEL/TEST**.25 
IF(tEST-.05)  60,  60,  70 
DEL=MIN1F  (  Wl  (211),  DEL/ TEST--- .  25  ) 

SELECTION  OF  STEP.  SIZE  FROM  ERROR  TEST 
DO  7  5  1=1,  2  40 
Dl  (  I  )=D3(  I  ) 
Wl  (  I  )  =W3  (  I  )• 

STORAGE  OF  INTEGRATED  VALUES  FOR  NEXT  INTEGRATION  STEP 
Wl  (223)  =  TEST--.  000003 
I F ( Wl ( 175 ) -FLAG )  80,  90,  80 
Wl (212) =W1 ( 175) -FLAG 
W  1  (  1  7  5  )  =  F  L  A  G 

STORAGE  OF  THRUST  LIMIT  SWITCH  POINT 
CALL  D  E  R  V ( W 1 ,  D 1 ) 
'■M  (2  12  )  =0.0 
RETURN 
END 
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1" 

]  5 


2  0 


RUN 

LAE 
L  T^ 
SUB 


d  i  y. 


<UT 
EL 
T 

ROUTINE 
SUBRO 
LET  X 
X 
(  XO  ) 
(XI) 
(  X2  ) 
(  X3  ) 
DEL 
ENS  I  ON 
)  =  .5 


DEL) 
INTEGRAT ION 


RUNKUT(W1»  Dl 
UTINE  PERFORMS 
*  =  DX/DT 
X»T  J 

XO,  TO) 
( XO  +  .5 (DEL 
(XO  +  .5(DEL 
XO  +  ( DEL  T ) ( X2 ) 
EL  T  )  (  (XO)*  +  2 


USING  KUTTAS  SIMPSONS  RULE 


*  =  F( 

*  =  F  ( 
-«-  =  F(   (XO  +  .5  (DEL  T)(XO)*)»  (TO  + 

*  =  F(  (XO  -r  .5  (DEL  T)(Xl)*)j  (TO  + 

*  =  F ( X  0  +  (DEL  T  )  (  X  2  )  *  )  >  (TO  +  DEL 
X  =  (DEL  T)(  (XO)*  +  2 (XI  J*  +  2(X2)* 
WI(240)>  D1I163),  W(240)»  D(168»  4)» 


.5  (  DEL 
.5 (DEL 
T)   ) 

+  (X3) 

CO) 


T)  ) 
T)  ) 

*  )  /  6 


C(2  J 

C  (  3  ) 

DO  5 

W  (  I  -r 

D  (  I  , 
DO  1 
F  =  DE 
TO  1 
'•'(  I  ) 
CALL 
E  =  DE 
DO  2 
W  1  (  I 
CALL 
RFTU 
END 


=  .5 
=  1.0 

1  =  1. 
7  2  )  =  W 

1  )=D 

5  L  =  l 
L*C(L 

0  1  =  ] 
=  W1  (  I 

DERV 
L  /  6  .  0 
0  I  -  1 
)=W1  ( 

DERV 
RN 


168 

i  (  1+72 ) 

KI) 

»  3 

) 

,  168 

) +E*D( I 

(  W  ,  D  (  1 


L) 
L  +  l  ) 


,  168 

I  )+E*(D( I » 

(  W 1  ,  D 1  ) 


1)+2.0*<D(I,  2)+D(I,  3))+D(I»  4)) 
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5 

10 


11 


12 


14 


15 


D   20 


25 


35 


INVERT 

LABEL 

LIST 

SUBROUTINE  INVERT(N,  QO ) 

SUBROUTINE  INVERTS  MATRIX  BY  SIMULTANEOUS  DOUBLE  PRECISION 

ROW  REDUCTION  OF  THE  MATRIX  TO  IDEM  AND  IDEM  TO  THE  INVERSE 
DIMENSION  00(10,  10),  Q(10,  20) 
DO  10  1=1,  10 
DO  5  J=l ,  10 

J) =QQ( I  ,  J  ) 

J+10  )  =0.0 

J  +  2  0  )  =0.0 

J+30 ) =0.0 
0 


I ) )-ABSF(Q( J,  I ) ) )  11,  14, 
LARGEST  ELEMENT  IN  COLUMN 


K+10) 


Q  (  I 

0(  I 

0(  I 

0(  I 

0(1,  1  +  10  )  =  1. 

DO  30  1=1,  N 

DO  14  J=I ,  N 

IF(ABSF(Q( I , 

TEST  FOR 
DO  12  K=l,  N 
S=0(J,  K) 
Q(J,  K }  =  Q  (  I  ,  K) 
0(1,  K  )  =  S 
S=Q(J,  K+10) 
0 (J,  K+10  )  =0(  I  , 
0(1,  K+10) =S 

TRANSFER  ROW 
CONTINUE 
D I V  =  Q(  I  ,  I  ) 
DO  15  J=l,  N 
0(1,  J ) =Q( I ,  J ) /DIV 
0(1,  J+10 ) =Q(  I  ,  J+10 )/DlV 

DIVISION  BY  DIAGONAL 
DO  30  J=l,  N 
IF(I-J)  20,  30,  2  0 
DI V=Q( J,  I  ) 
DO  2  5  K=l,  N 

Q(J,  K)=0(J,  K)-Q(I,  K )*DI V 
0(J,  K+10)=0(J,  K+10)-Q(I,  K+10) 

DIAGONALIZATION  OF  MATRIX 
CONTINUE 
DO  35  1=1,  N 
DO  35  J=l,  N 
00 ( I ,  J)=Q( I ,  J+10 ) 
RETURN 
END 


14 


OF  LARGEST  ELEMENT  TO  FIRST  ROW 


ELEMENTS 


•DIV 


45 


TOTAL 


239* 
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Dl  ) 


3  )*>.'(  2  ,  6  )  ) 


0(5, 

0(5, 

n  ( <s , 

D(5> 

D(5» 

0(5, 

c 

W  (  6  » 

c 

W  (  6  » 

W  (  7  , 

c 

DO    1 

'/J  (  6  , 

^_ 

D(l, 

D(  1» 

c 

D(2t 

10    D(2, 

C 

D(l. 

C 

D(2  » 

c 

0(5, 

<" 
L 

00    2 

0(3, 

c 

.-  (  4  , 

1.0) 


DERV 

LABEL 

LIST 

SUBROUTINE  DFRVt 

LINEAR  OPTIMUM 

SUBROUTINE  ESTABLISHES  DIFFERENTIAL  EQUATIONS  OF  SYSTEM 
DIMENSION  W(12»  20),  D(12»  14),  DK168) 
R  =  SGRTF(W(1,  1)*#2+W(1»  2)**2+W(ltv  3)*#2) 
RH3=W(7»  17)./R**3 
RM5=3.0*RM3/R#*2 
RL=RM5*(W(1.  l)*W(2i  4)+W(l,  2)*W(2»  5 ) +W ( 1 , 

COMPUTATION  OF  GRAVITATIONAL  GRADIENT 
A2  =  SQRTF(  W(2  »  4)'**2+W(2»  5)**2+W(2»  6)**2) 

1  )  =1.0 

2)=W(6»  17)/W(1,  7)*W(7»  15) 

3)=D(5»  2)**2 

4)  =.5*A2-"-D(  5,  2) 

5) =0(5,  2)*W(1»  7) 

6)=D(5,  2)/W(li  7)*W(6»  17) 

DERIVATIVES  OF  PARAMETRIC  VARIABLES 

15  )  =R 

STORAGE  OF  R  FOR  ERROR  TEST 
16)=MAX1-F(A2*W(1»  7)/W(6,  17)*W(8»  15), 

16  )  =A2*W(  1  ,  7)#W(8,  15) 
COMPUTATION  OF  THRUST  SWITCHING  FUNCTION 

C  1=1,  3 

I  +  1  7  )  =  W  (  2  ,  I  +  3  )  *  D  (  5  ,  2  )  /  A  2 

COMPUTATION  OF  ACCELERATION  VECTOR 

I)=W(li  1+3) 

I+3)=W(6»  I+17)-RM3*W(1»  I) 

EQUATIONS  OF  MOTION 

I ) =RM3*W( 2,  1+3 )-RL*W( 1 ,  I ) 

1+3) =  -W( 2,  I  ) 

EULER  EQUATIONS  (MASS  INDEPENDENT) 

7)=-yJ(6,  17)/W(8,  16)*W(7»  15) 

MASS  RATE  EQUATION 

7  )=/,<(  2,  7)---D(5,  2)/W(8,  15) 

EULER  EQUATION  FOR  MASS 

7)=W(5»  7)*D(5,  2)/W(8»  16) 

ADJOINT  EQUATION  FOR  MASS  SENSITIVITY 
C  1=1,  6 

I  )=','.'(  1+5  ,  4)*W(6»  18)+W(I+5»  5J*W(6»  l9)+W(I+5» 

EQUATION  FOR  ETA  VECTOR 

I)=WU+5i  4)*W(2»  4)+W(I+5,  5)*W(2»  5)+W(I  +  5,  6 )  #W  (  2 » 

QUANTITY  USED  IN  M  STAR  MATRIX 
RL=RM5* ( W( I+5»  4)-W(l,  i)+W(I+5,  5 ) *W ( 1 >  2 
0(1+5,  7 ) =D( 3,  I ) /W( 1 ,  7) 

ADJOINT  EQUATIONS  FOR  MASS  DEPENDENCE 
00  2  0  J=l,  3 

0(1+5,  J)=RM3*W(  1+5,  J  +  3  )-RL-"-W(  1  ,  J) 
D (1+5,  J+3 )=-W( 1+5,  J ) 

ADJOINT  EQUATIONS  (MASS  INDEPENDENT) 


TOTAL 


6) 


(6,  20) 


6  ) 


+  W(  1  +  5,  6  )#W (  1 »  3 ) 


52 


150 


25 


R  L  =  (  W  (  2  »  1  )  #  W  (  2  ♦  4 )  +  W  (  2  ♦  2  )  *  W  (  2  »  5  )  +  W  (  2 

QUANTITY  USED  IN  M  STAR  MATRIX 


3)*W( 2*  6) )#A2*W( 1»  7 ) 


1  =  1 
U=I 


DO  30 

DO  3  0 

F  =  0  .  0 

DO  25  K=l, 

F  =  F+W(  1+5  » 


6 
6 

3 
K+3 


*W(J+5»  K+3) 


D ( 1+5 ,  J+7)=(F-D(3»  I)*D(3t  J)/D(5»  3))*D(5.  2J/A2 


3J 


DERIVATIVE  OF 
W (  1+5 »  J  +  7  )=W(  1+5  ♦ 

CORRECTION  TO 
W (U  +  5 ,  1+7 )=W(  1  +  5  > 
D( J+5,  1+7  )=D(  I+5» 


M  STAR  MATRIX 

J+7)+W(4»  I)#W(4»  J)/RL*W(6»  17 

M  STAR  FOR  VARIABLE  INTEGRATION 

J  +  7) 

J  +  7) 


*W(8 
TIME 


18) 


M  STAR  IS 
DO  40  1=1,  168 
D 1  (  I  )  =  D  (  I  »  1  ) 

RETURN 
END 


SYMMETRIC 


TOTAL 


18 
18* 


151 


EQUATIONS  OF  SYSTEM 


"  ( <s , 

W  (  7  » 

c 

00    1 

W  (  6  » 

r 

Dd, 

D(  1  . 

c 

0(2, 

10     0(2% 

c 

0(1, 

c 

0(2, 

c 

0(5, 

c 

OFRV 

LABEL  . 

LIST 

SUBROUTINE  DERV( W»  01 ) 

QUADRATIC  OPTIMUM 

SUBROUTINE  ESTABLISHES  DIFFERENTIAL 
DIMENSION  W(12t  20),  D(12,  14),  DK163) 
W  (  «  ,  1  5  )  =  1  .  0 

R=SQRTF(W(1,  1  )**2+W( 1 ,  2)**2+W<l»  3)**2) 
RM3=W(7,  17)/R**3 
RM5=3.0*RM3/R**2 
RL  =  RM5-(  .'.'(  1  ,  1)*W(2»  4)+W(l,  2)#W(<2,  5  ) +W  (  1  ,  3)*W(2»  6)) 

COMPUTATION  OF  GRAVITATIONAL  GRADIENT 
A2=SQRTF( W(2 »  4)**2+W(2»  5)*#2+W(2,  6)**2) 
GMASS  =  W ( 2 ,  7 ) *W (  1  ,  7 ) *#2/W ( 8  >  17) 
G AM=MAX 1 F ( A2/W ( 6  »  17 ) *W ( 1 »  7 ) -GMASS ,  1.0) 

THRUST  LIMIT  SWITCHING  FUNCTION 
STAR=W(7»  15)#(W(1»  7)/W(6,  17))**2 

QUANTITY  USED  IN  M  STAR  MATRIX 
0(5,  1  ) =  1  .  C 

0(5,  2) =A2/( GAM+GMASS) 
0(5,  3 ) = D ( 5 ,  2 ) * *  2 
0(5,  4) =.5*GAM*D( 5 ,  3) 
0(5,  5 )  =  0 ( 5  ,  2 } * W ( 1 ,  7 ) 
0(5,  -S)=0(5,  2)/W(l,  7)*W(6,  17) 

DERIVATIVES  OF  PARAMETRIC  VARIABLES 
if.' (6»  15)=R 

STORAGE  OF  R  FOR  ERROR  TEST 

16 ) =GAM 

1 6 ) =  A2*W (  1  »  7 ) -GMASS*W ( 6 ,  17) 

STORAGE  OF  THRUST  LIMIT  FUNCTION 
1  =  1  ,  3 

I  +  1  7  )='•■:(  2  -.  1+3  )/ (GAM  +  GMASS) 

COMPUTATION  OF  ACCELERATION  VECTOR 

I  ) =W(  1,  1+3  ) 

I+3)=W(6»  I+17)-RM3*W( 1 ,  I) 

EQUATIONS  OF  MOTION 

I  )  =RM3*W< 2  ,  1+3 )-RL*W( 1 ,  I  ) 

1+3 ) =-W( 2  ,  I  ) 

EULER  EQUATIONS  (MASS  INDEPENDENT) 

7)=-D(5,  3)*W(1,  7)**2/W(8»  171/2.0 

MASS  RATE  EQUATION 

7)=W(1,  7)*D(5,  3)/W(8»  17)*W(2»  7) 

EULER  -EQUATION  FOR  MASS 

7  )  =W  (  1 ,  7  )  *D ( 5 ,  3  )  /';.'(  6  ,  17  )  #W  (  5  ,  7  ) 

ADJOINT  EQUATION  FOR  MASS  SENSITIVITY 

oo  z:   i=i,  6 

0(3,  [)=W(I+5,  4)---W(6,  18)+W(I+5,  5)*W(6»  19  1+WU+5,  6 )  *W  ( 

EQUATION  FOR  ETA  VECTOR 
RL=RM5*(W(  1+5,  4  )  - W  (  1  ,  l)+W(I+5,  5  )  *W  (  1  ,  2J+WII  +  5,  6)*Wtl» 
0(1+5,  7)=D(3,  I  )/W(  1  ,  7  ) 

ADJOINT  EQUATIONS  OF  MASS  DEPENDENCE 


i,  20) 
3  )  ) 


TOTAL 


52- 
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20 


25 


DO  20  J=l»  3 

D  (  1+5,  J) =RM3#W( I+5»  J- 

0(1+5,  J  +  3  )  =  -  'a'  (  I  +  5  ,  J  ) 
ADJOINT  EQUATIONS 


3)-RL*W(l,  J) 


(MASS  INDEPENDENT) 


DO 
DO 


1  =  1 
J=I 


3  0 

30 
F=G.O 

DC  25  K=l, 
F=F+W( 1+5 , 
D  (  1+5,  J  +  7 


3 

,<  +  3)*W(J  +  5»  K  +  3) 

=(F-D(3,  I)*D(3»  J)*STAR)/GAM 

DERIVATIVE  OF  M  STAR  MATRIX 
D( J+5,  I+7)=D( 1+5,  J+7) 

M  STAR  IS  SYMMETRIC 
DO  40  1=1 ,  168 
Dl  (  I  )=D(  I  ,  1 ) 
RETURN 
END 


TOTAL 


17 
17* 
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APPENDIX  H 
NUMERICAL  RESULTS 

H.  1  Introduction 

In  this  appendix  the  computer  output  data  from  a  sample  run  of 
each  control  mode  are  reproduced.     The  data  sheets  are  in  general 
self  explanatory  except  perhaps  for  the  format  of  the   A    matrix  and 
the  units. 

All  parameters  and  variables  derived  in  the  linear  theory  are 
printed  out.     In  addition,  the  important  variables  are  plotted  in  Figures 
H-a  through  H-r. 

H.  2  Format  and  Units  of  the  Data 

The  first  pages  of  each  example  contain: 


1. 
2. 
3. 

4. 
5. 
6. 

7. 


The  time  in  days  from  initial  point, 


-4 


The  acceleration  program  in  units  of  10      g  , 

-4 
The  acceleration  magnitude  in  units  of  10      g  , 

The  position  in  A.  u.  , 

The  velocity  in  A.  u.     per  day, 

The  normalized  mass, 

The  switching  parameter,    y  . 


The  coordinates  are  computational  coordinates  in  all  cases. 

Following  the  state  and  acceleration  data  are  the  terminal  values 
of  interest  printed  as  column  vectors.     The  units  are  A.  u.    and  A.  u. 
per  day. 

The  elements  of  the  M*  matrix,  which  are  the  next  set  of  numbers, 
are  not  particularly  interesting  and  may  be  disregarded. 

Adjacent  to  the  M*  matrix  is    the    cost    quantity  J  with  units  of 
(A.  u.  )     per  (day)    . 
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The  next  set  of  pages  contain  the  elements  of  the    A  matrix  for 
both  optimal  acceleration  and  optimal  thrust  programs.     For  each  time 
point,  the  array  consists  of  50  elements  which  are  interpreted  as  par- 
tial derivatives  or  influence  coefficients.     The  array  is  ordered  as  in 
the  equation 

6sf  =    A(t)  6st 

The  first  49  elements  are  the  adjoint  set  for  an  optimal  thrust  program. 
The  adjoint  set  for  an  optimal  acceleration  program  is  obtained  by  re- 
placing the  seventh  column  with  zeros  for  the  first  six  elements,  and 
replacing  the  1  with  the  element  below  it.     That  is:     the  array 


All       A12 


A 
O 


21 
T 


A. 
O 


22 

T 


Al3 
A23 


A. 


33 


yields 


A  (optimal  thrust) 


A 
A 

O 


11 

21 
T 


A12     A 13 

A22     A23 


o 


and 


A    (optimal  acceleration) 


An 
A 


o 


21 
T 


A12 
A22 

oT 


o 
o 


A 


33 


The  applicable  submatrices,  interpreted  as  partial  derivatives,  are 
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dL{ 

a£f 

dr. 

8v. 

—  t 

— t 

*Zl 

8^f 

tt 

^t 

3m 

3mf 

8r. 

8i't 

— t 

9r 


3m, 


3v 


3m, 


3m, 


3m, 


The  units  are  A.  u.    and  A.  u.    per  day  for  position  and  velocity.     The 
seventh  column  is  with  respect  to  a  100%  change  in  mass.     For  example 


3r 
the  units  of      —  f     are 

3m, 


A.  u. 


100%  change  in  mass 


The  first  set  of  data  is  for  unrestricted  VSI  control.     The  second 

set  is  for  CSI  control  with  a     =  1.  2  X  10~4g   . 

o  &o 
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2.2 

2.  I 
2.0 

1.9 
1-8 
1.7  - 
1.6  - 


0-5 
0-4 

0-3 

0-2  ~ 


VARIABLE    SPECIFIC  IMPULSE 


CONSTANT 
SPECIFIC 
IMPULSE 


0        10       20      30      40      50      60      70       80      90      100      110      120      130     140     150 

TIME  (DAYS) 

Fig.  H-a.    Acceleration  schedule  for  150-day  Earth-Mars  transfer. 
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Fig.  H-d.    Sensitivity  of  rx  (tf)  to  position  variations. 
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Fig.  H-e.    Sensitivity  of  rx  (tf)  to  velocity  variations. 
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23  3r, 
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Fig.  H-f.    Sensitivity  of  ry  (tf)  to  position  variations. 
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C  D 

A24  = 

ary  (tf) 

(METERS) 

24 

*vx 

(METERS/SECOND) 

22 

—       \ 

K*- 

ary(tf) 

0Vy 

20 

\                  A26  = 

ary(tf) 
dv2 

1-8 

1-6 

1.4 

1.2 

\  A25  x  I0"7 

10 

0.8 

0.6 

V\A  26  XKT5 

0.4 

0.2 

0 

A24  X  I0~7 

-0.2 
-n  4 

1 

1             1             1 

i        I        1 

1          1           1          1 

1           1 

0         10        20       30      40       50      60        70      80       90      100      110       120      130      140     150 

TIME  (DAYS) 

Fig.  H-g.    Sensitivity  of  ry  (tf)  to  velocity  variations. 
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Fig.  H-h.    Sensitivity  of  rz  (tf)  to  position  variations. 
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Fig.  H-i.    Sensitivity  of  rz  (tf)  to  velocity  variations. 
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Fig.  H-j.    Sensitivity  of  vx  (tf)  to  position  variations. 
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Fig.  H-k.    Sensitivity  of  vx  (t{)  to  velocity  variations. 
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Fig.  H-l.    Sensitivity  of  vy  (tf)  to  position  variations. 
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Fig.  H-m.    Sensitivity  of  v    (tf)  to  velocity  variations. 
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Fig.  H-n.    Sensitivity  of  vz  (tf)  to  position  variations. 
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Fig.  H-o.    Sensitivity  of  vz  (tf)  to  velocity  variations. 
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Fig.  H-p.    Sensitivity  of  r  (tf)  to  mass  variation  for  optimal  thrust  program. 
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Fig.  H-q.    Sensitivity  of  v  (tf)  to  mass  variation  for  optimal  thrust  program. 
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Fig.  H-r.    Sensitivity  of  m(tf)  to  mass  variations  for  optimal  acceleration  program. 
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